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ON ■ In this review the problem of statistical description of isolated quantum systems of interacting 

particles is discussed. Main attention is paid to a recently developed approach which is based 
on chaotic properties of compound states in the basis of non-interacting particles. In order to 
demonstrate the most important aspects of this approach, the matrix model of two-body random 
interaction between Fermi-particles has been used. Different problems have been considered such 
as the onset of chaos and statistical equilibrium, the relation between the structure of eigenstates 
and distribution of occupation numbers, the emergence of the Fermi-Dirac distribution in isolated 
systems of finite number of particles and many others. The application of the approach to dynamical 

l ' systems with the classical limit is discussed as well. 
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I. INTRODUCTION 

> \ 

, Until recently, the quantum chaos theory was mainly related to few-body physics. On the other hand, in real 
physical systems such as many-electron atoms and heavy nuclei, the origin of complex behavior is quite strong 
interaction between many particles. To deal with such systems, famous statistical approach has been developed which 
is based on the Random Matrix Theory (RMT) (see, for example, The main idea of this approach is to forget 

about a detailed description of the motion and to treat these systems statistically having in mind that the interaction 
between particles is so complex and strong that generic properties are expected to emerge. Simplest models of the 
RMT are full random matrices of a given symmetry, the latter was shown to have a direct link with underlying 
symmetries of physical systems. 

One of the main results of the RMT is the prediction of a specific kind of correlations in the energy spectra of 
complex quantum systems. Among many characteristics of these correlations, the most popular is the distribution 
j-^, of spacings between nearest energy levels in the spectra. Exact analytical expression of this distribution is very 
' complicated, instead, one uses the so-called Wigner-Dyson (WD) surmise (quite simple expression which gives a very 
good approximation to the exact result). A distinctive property of this WD-distribution is the repulsion between 
neighboring levels in the spectra, the degree of this repulsion (linear, quadratic or quartic) depends on the symmetry 
of random matrices. In fact, such type of repulsion was observed experimentally very long ago (first experimental 
observation is reported in Ref. || for the energy spectra of heavy nuclei), and discussed in many theoretical works. 

After this prediction of the RMT, the WD-distribution has been confirmed to occur in heavy nuclei and many- 
electron atoms, see references in |^||. Later on, it was found also in dynamical systems with chaotic behavior in 
the classical limit, famous examples are the so-called billiards (see for example, Q). As a result, it was understood 
that chaotic properties of quantum systems are generic for both disordered models (when the randomness of matrix 
elements is postulated from the beginning) and dynamical systems for which the pseudo-randomness appears as a 
result of special conditions, the latter are convenient to explain by comparing with the classical limit. Thus, one 
can say that limiting properties of quantum chaos in the case when all regular dynamical effects are neglected, are 
described by the RMT. 

As one can see, the RMT can give a proper description of a system (mainly, the properties of energy spectra) 
only locally, in a restricted region of energy spectra. Indeed, the RMT can not give any global energy dependence 
neither for the energy spectrum nor for eigenstates, it is parameter-independent theory. In this sense, the conventional 
RMT (ensembles of fully random matrices) is very restricted and, for example, it can not be directly applied to such 
important phenomena as the localization of eigenstates in disordered models. This is why new approaches in the 
RMT have been developed by imposing internal structure of random matrices. The most known example is the 
so-called Band Random Matrices (BRM), or random matrices with a band- like structure (see, for example, R], l| and 
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references therein). Inside the band, the matrix elements are assumed to be random and independent, and outside 
the band matrix elements are set to zero. The BRM-ensemble has been recently studied in details both numerically 
and theoretically, and much is now known about the structure of eigenstates and spectrum statistics for both infinite 
and finite matrices. The application of this kind of matrices in physics is very broad. In particular, they have been 
used to describe dynamical localization in dynamical systems with time-periodic perturbation (paradigmatic model 
is the Kicked Rotor @]), and the localization of eigenstates in quasi-ld disordered models in solid state physics |§. 

Another class of band random matrices has been introduced very long ago by Wigner in Ref. ||. The structure 
of these Wigner Band Random Matrices (WBRM) is characterized by the leading diagonal with reordered (in an 
increasing way) values, plus random and independent off-diagonal elements inside the band of size b . This type of 
matrices is much closer to physical realistic systems, compared to the standard full random matrices. One can suggest 
that the original motivation of Wigner for the study of these matrices was a close correspondence to Hamiltonians 
of complex nuclei, which are typically described by the mean-field part Hq (leading diagonal) and the residual 
interaction V of the finite energy range (off-diagonal matrix elements inside the band). The main interest of Wigner 
was the quantity which nowadays is known as the strength function or local density of states (LDOS). This quantity is 
extremely important when describing the spread of energy, initially concentrated in a specific state of the unperturbed 
Hamiltonian Hq , between all other states due to the internal interaction V . In particular, it was analytically shown 
that this function has the form of the Lorentzian if interaction is sufficiently (but not extremely) large. This result is 
fundamental and it is used in very different applications. 

However, in spite of the success of the standard RMT and its modern developments, the strong assumption of 
randomness of matrix elements, as well as a specific (band) structure of matrices do not allow to relate such matrices 
directly to realistic many-particle Hamiltonians. One of the important reasons is that the underlying structure of 
realistic Hamiltonian matrices results from the single-particle spectrum and two-body interaction between particles. 
In order to understand the role of k— body (random) interaction, a new ensemble of matrices has been suggested 
(see, for example, |l(],[ll],|3| and references therein ). In this approach the matrices arise as a result of construction 
from the single-particle basis, provided random character of the k— body interaction. The study of this ensemble of 
matrices have shown that even when the interaction V is very strong and the influence of the leading diagonal may be 
neglected, there are serious differences from the standard RMT. In particular, it was discovered that for the two-body 
interaction, k — 2 , the spectral fluctuations are different from those predicted by the RMT, although the distribution 
of spacing between nearest levels has the form similar to the WD-distribution. It was shown that full random matrices 
occur when the rank of interaction is very large, k — > oo . 

Due to very serious mathematical problems these two-body random interaction (TBRI) matrices were forgotten for 
quite a long time, and only recently they have been used in the context of quantum chaos. In these lectures, the author- 
gives a review of recent results obtained for the TBRI-model in collaboration with the co-authors. Main attention 
is paid to a novel approach developed in 12-lq|, which is based on the chaotic structure of eigenstates in a given 
basis of unperturbed many-particle states. This approach allows to relate statistical properties of exact eigenstates 
in many-body representation directly to properties of single-particle operators, in the first line, to the occupation 
number distribution of single-particle states. 

The structure of the paper is as follows. In the next Section 2.1 the structure of the TBRI-matrices is discussed 
in details. It was explained how these matrices are constructed from single-particle states and what are properties of 
matrix elements in many-body representation, also, the comparison is made with full random matrices. A particular 
point is that in spite of complete randomness and independence of two-body matrix elements, off-diagonal matrix 
elements of many-particle Hamiltonian matrices have underlying correlations which are due to a two-body nature of 
interaction. 

In Section 2.2 generic properties of density of states and level spacing distribution are briefly discussed, although 
they are not the main interest of the study. Section 2.3 deals with the structure of exact eigenstates in dependence 
on the interaction strength and excitation energy. The notion of the average shape of the eigenstates (F— function) 
is introduced and discussed in details since the structure of chaotic eigenstates plays a basic role in the approach. 
In next Section 2.4 the structure of the strength function is considered and compared with that of exact eigenstates. 
Main attention is paid to the conditions under which this form is the Lorentzian, and how this form changes with 
an increase of interaction strength. In Section 2.5 very recent analytical results are discussed, which arc obtained 
for the shape of the strength function in the TBRI-model for any strength of interaction. Section 2.6 is devoted 
to non-statistical properties of this model. Specifically, it is shown that some quantities can not be described in a 
statistical way, in spite of completely random character of the two-body interaction. 

Next Section 3.1 starts with the discussion of the basic relation between the structure of exact eigenstates and the 
distribution of occupation numbers (DON). It is shown that the average shape of eigenstates plays a crucial role for 
this distribution and practically determines mean values of single-particle operators. In next Section 3.2 the relevance 
of the DON for isolated systems, to the standard canonical distribution is discussed. It is shown how the canonical 
distribution emerges in isolated systems with an increase of the number of interacting particles. In next Section 
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3.3 the problem of the Fermi-Dirac distribution is considered for the TBRI-model of interacting Fermi-particles, in 
particular, conditions under which this distribution occurs in isolated systems are analyzed. An important problem of 
the statistical description of the occupation number distribution is considered in Section 3.4. An analytical approach 
has been developed in order to obtain the DON in the case of statistical equilibrium which results from the chaotic 
structure of eigenstates. This approach is valid even for small number of particles, in the case when the DON differs 
from the Fermi-Dirac distribution. In next Section 3.5 another approach is suggested for the description of the DON, 
in the case when its form is of the Fermi-Dirac type. It was shown how the DON can be obtained by a proper 
renormalization of the total energy of a system, originated from the interaction between particles. General discussion 
of the meaning of temperature in isolated systems of finite number of particles is the content of Section 3.6. The main 
point is that for small number of particles different definitions of temperature give different results. Therefore, is it of 
interest to compare these definitions and to understand their meanings, if any. Finally, in Section 3.7 main results are 
summarized for the transition to chaos and equilibrium in the TBRI-model in dependence on the interaction strength. 
In last Section 4 it is briefly shown how the developed approach can be applied to dynamical systems. 



II. TWO-BODY RANDOM INTERACTION MODEL 



A. Description of the model 

1. Many-body Hamiltonian 

The model we discuss here deals with Hamiltonians which can be separated in two parts, 

H = H a + V (1) 

where Ho describes the "unperturbed " part and V stands for the interaction between particles or between different 
degrees of freedom. In order to study statistical properties of such models we assume in the following that the 
interaction is completely random. In contrast with standard approach of the RMT where matrix elements of V are 
taken as random variables, we would like to keep an important physical property of real systems and to take into 
account that the interaction is of the two-body nature. Therefore, we start with the single-particle Hamiltonian which 
refers n non-interacting particles occupying m single-particles levels, and assume that the matrix elements of the 
two-body interaction V SlS2S3Si are independent random variables. Here, the indices si, S2, S3, S4 indicate initial (si, S3) 
and final (s 2 , s 4 ) single-particle states coupled by the interaction. 

In what follows we consider Fermi-particles, however, the approach can be easily extended to Bose-particles |]l7| , |l8| . 
Therefore, in the Slater determinant basis the unperturbed part has simple form 

H = ^2 e-s a\a s (2) 

and the perturbation can be represented as 

v = \^2Vs lS2 s 3 s i al i al 2 a S3 a Si . (3) 

Here e s is the energy of a particle, corresponding to the single-particle state |s) and a\. , a s . are creation-annihilation 
operators. With these notations, exact eigenstates \i) of the total Hamiltonian H {compound states) can be expressed 
in terms of eigenstates k of the unperturbed part Ho (basis states) as follows 

|i>=E C * ) l fe )> \k)=al 1 ...al\0) (4) 

k 

where is the k — th component of the compound state \i) in the unperturbed basis. These components determine 
the important quantity which will be discussed in great details below, the occupation numbers n s , 

n s = (i\n s \i) = Y,\C®\ 2 (k\n s \k) (5) 

with n s — a|a s as the occupation number operator. For Fermi-particles the occupation number ni = (k\n s \k) is 
equal to 1 or depending on whether any of the particles in the basis state | fc) occupies or not the single-particle state 
Is) . 
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As one can see, our model is described in terms of many-particles basis states |fc) and exact (compound) states \i) 
which are constructed from the single-particle states \s) and two-body matrix elements V S1S2S3S4 . In what follows we 
assume that the basis states are reordered in an increasing way for the total energy = 2_, s e s n s with an increase 
of the index k = 1, . . . , N . This way of the ordering of the unperturbed basis is crucial for the analytical description 
of chaotic compound states which are formed by the interaction between many basis states, see below. The size N of 
the basis for the many-particles Hamiltonian H can be found from the combinatorics: if any single-particle state can 
be occupied by one particle only, one can get 



where to is the number of single-particle states (orbitals). The latter estimate in (^|) shows that total number of 
many-particle states increases very fast (exponentially) with an increase of the number of particles and orbitals. For 
example, for m = 11 and n = 4 the size of the H— matrix is N = 330 . The model with these parameters has been 
studied in great details in fl^-|l5f| and compared with direct computations of the Ce atom fl9|-pl"|| . For this atom there 
are n = 4 valence electrons, and the core can be effectively described by the Hartree-Fock method. This method has 
been used []l9|-f2l| in order to calculate the basis set of single-particle relativistic states with energies e s as well as the 
matrix elements Vs lS2 s 3 s i of interaction between valence electrons. The Ce atom is known to have good statistical 
properties, and this was the reason in ||l2]-|l5| to compare direct calculations with the simplest model (Q,^|h we are 
going to discuss. In spite of the fact that this model does not take into account the momentum (it depends only on the 
energy and for this reason can be treated as "zero-dimensional " , it turns out to be instructive for the comparison of 
the real (dynamical) Ce atom with the two-body random interaction model (TBRI-model) described above (for details 
see [J12H15|). In what follows, for the single-particle energies e s we take a non-degenerate spectrum with constant 
mean level spacing do = < e s+ \ — e s > which, without the loss of generality can be taken do = 1 . The unperturbed 
single-particle spectrum has been chosen at random, or according to the expression e s = do(l + 1/s) (the results 
are statistically the same). As one can see, the model is defined by four parameters, to, n, do and V 2 = \V s ^ S2S3S4 j 
which is the variance of two-body random matrix elements (we assume that the distribution of these elements is the 
Gaussian with the zero mean). 

The Hamiltonian (Q) with (||) and (^) is of general form, it appears in many physical applications such as when 
describing complex atoms, nuclei, atomic clusters etc. In fact, the form of H discussed above is known as the mean 
field approximation for complex quantum systems of interacting particles. In this description, the unperturbed part 
Hq represents the zero-order mean field for the excited states with the ground state E\ , and the residual two-body 
interaction is given by V . Therefore, the single-particle levels e s in such applications are, in fact, renormalized 
quasi-particle energies (see details, for example, in (22)). The considered here model does not take into account such 
physical effects as momentum dependence, pairing effects and others, however, it contains the main effects of quantum 
chaos and is very effective for the understanding generic features of complex systems. Moreover, as was pointed out, 
the approach we discuss below, can be extended for dynamical quantum systems which exhibits complex behavior. 

From the view point of dynamical systems when the complexity of the behavior appears as a result of dynamical 
chaos (both for the systems with or without the classical limit) , the separation of the total Hamiltonian in two parts is 
well defined physical procedure. Specifically, the approach is expected to be valid if the second part V is as "random 
" as possible. In other words, one should find such separation that the perturbation has no strong regular part, the 
latter should be embedded into the "unperturbed " Hamiltonian Hq . In this way, the compound states may be 
treated, for sufficiently strong perturbation, as chaotic superposition of simple basis states. As is well known, such 
a situation is typical for many-electron atoms and heavy nuclei. Indeed, the number of basis states (number N pc of 
principal components ) is known to be about 10 — 10 6 for excited nuclei, and ~ 100 in excited rare-earth or actinide 
atoms. 

By assuming the complete randomness of two-body elements in our models, we avoid the influence of any regular 
effects. Therefore, our goal is to explore statistical properties of the model in its strongest "chaotic " limit. One 
should stress that the answer is far from being trivial since the many-body Hamiltonian H turns out to be quite 
different from standard random matrix ensembles of the RMT and, as will be shown, some important quantities can 
not be described statistically. 

It should be noted that the TBRI-model we discuss here, for the first time was analyzed long ago (see |Io| , |ri| . 
also, the review || and references therein). The original interest was related to the fact that for standard random 
matrices the density of states has the famous semicircle form, in contrast with physical systems for which the RMT 
was addressed (complex nuclei). Moreover, the density of states for many-body systems increases very fast, and 
has nothing to do with the semicircle even on a local scale. Therefore, the natural question is to understand what 
is missed in the RMT. For this reason, another ensemble of matrices has been suggested which takes into account 
fc— body interaction between particles. The theoretical study has been shown that, indeed, fully random matrices of 
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the kind considered in the standard RMT, formally correspond to the case k — > oo . On the other hand, with decrease 
of k , the density of states tends to the Gaussian form. The latter form is closer to reality and on a local scale, the 
density may be treated as the realistic one. Another question which has been under close investigation, is the role 
of fc-body interaction on the spectrum statistics, see discussion below. For some reason, these studies have not been 
extended until recently, when the role of two-body interaction in different applications was questioned in the context 
of quantum chaos. Unlike the previous studies, below we pay the main attention to chaotic properties of eigenstates, 
and to the problem of how these properties can be linked to the properties of single-particle operators, such as the 
occupation number distribution. 



2. Structure of the Hamiltonian matrix 

Let us start with the structure of the Hamiltonian H in the chosen basis. For this, one should construct matrix 
elements Hij = (i\H\j) which correspond to the coupling between basis states \i) and \j) due to the interaction 
V. One can immediately see that the number N2 — m 2 (m — l) 2 /2of independent matrix elements V SlS2S3S4 of the 
two-body interaction is much less than the total number N(N + l)/2 of the (symmetric) matrix elements . It is 
very important that due to a two-body character of the interaction, the matrix elements are non-zero only when 
basis states \i) and \j) differ by no more than two occupied single-particle states. In order to count the total number 
K of non-zero matrix elements for the fixed i, we separately count the numbers K , K\ , K 2 of non-zero matrix 
elements which correspond to the transition between the basis states which differ by the positions of none, one and 
two particles, respectively, 

Kq = 1, K\ — n(m — n), K 2 = -n(n — l)(m — n)(m — n— 1). (7) 

As a result, the total number K in each line of the matrix is 

K = K + Ki + K 2 = 1 + n(m - n) + ~n(n - l)(m - n)(m - n - 1) m ^n 2 m 2 (8) 

where the last estimate is given for large number of particles and orbitals, 1 « n C m . Comparing K with ./V 
defined by Eq.(^) , one can see that the matrix H is sparse (only for n = 2 there is no forbidden transitions and the 
matrix is full). 



50 100 150 200 250 300 

n 2 

FIG. 1. Sparsity of the Hamiltonian matrix H ni: , l2 for n = 4 particles, m = 11 orbitals. Black points are non-zero matrix 
elements. 



The presence of many zeros in the matrix means that, in a sense, there are strong correlations between matrix 
elements since the position of zeros are fixed for any random choice of the two-body random matrix elements. It 
should by pointed out that the sparsity increases with an increase of the number of particles, therefore, the more 
particles (and orbitals), the less relative number of non-zero elements. This fact is important for the description of 
such systems by the random matrix approach. One should remind that in the standard RMT, the sparsity is not 
taken into account at all. 

The sparsity of the matrix for m = 11 and n = 4 is shown in Fig.l where black points correspond to non-zero 
elements. First, one can see that the density of zero elements increases when moving away from the principal diagonal. 
Second, the positions of non-zero matrix elements are correlated, there are some curves along which the density is 
high, this reflects the two-body nature of interaction. In contrast with full random matrices of the standard RMT, 
the influence of the off-diagonal elements depends on the distance from the principal diagonal. To illustrate this 
peculiarity, we averaged the modulo of the off-diagonal matrix elements over blocks of the fixed size 10 x 10 in such 
a way that instead of the matrix of the size N x N we have the reduced size N/10 x N/10 . The result is shown in 
Fig. 2 where only off-diagonal terms are presented. As one can see, the amplitude of these effective matrix elements 
decreases when moving away from the diagonal. This means that the effective intensity of the off-diagonal terms 
far from the diagonal is less than of those close to the diagonal. In some sense, one can treat the structure of the 
Hamiltonian as the band-like, although it is clear that the amplitude of the averaged matrix elements decays quite 
slowly. 




FIG. 2. Shape of the Hamiltonian matrix (without leading diagonal Hi t i). The matrix for n = 4, m — 11, Vq = 0.12, do = 1.0 
has been divided into blocks of size 10 x 10 and the sum |J? n , m | = y\ . \Hi,j\ has been computed inside each block (n, m) . 
This type of the average allows to see effective intensity of the off-diagonal matrix elements. 



3. Correlations in off-diagonal matrix elements 

Now, let us analyze non-zero off-diagonal elements of the matrix . According to the definition (0) , any of these 
elements is just a sum of one or more two-body matrix elements. When the basis states \i) and \j) differ by one 
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occupied orbital, \j) = a\ a si \i) , the matrix element Hij is the sum of n — 1 two-body matrix elements, 



n-2 



(9) 



Here the last equality is given in order to show that among other matrix elements Hy y there are such elements which 
differ from by the last term only and the sum of n — 2 is exactly the same fl3| , 



n—l n-2 



(10) 



with a ^ [3 . This happens for those basis states whose many-particle energies differ by the energy difference 
corresponding to the move one particle from the orbital a to the orbital f3 , 



E, 



Ei — Eji — Ej — ep — e a 



(11) 



One can see that these matrix elements, strictly speaking, can not be treated as completely independent variables. 
Therefore, if one averages over the ensemble of matrices which constructed from different sets of two-body random 
elements, the correlations for such elements remain, (HijHi/ji) ^ . 

The more striking result arises when considering matrix elements which correspond to the coupling between 



those basis states \i) and \j) which differ by two occupied orbitals, |j) 



These matrix elements are 



equal to the corresponding single two-body matrix elements, in other words, there is only one term in the sum in 
Eq. (^), Hij = V Sllll S2fJ _ 2 . These matrix elements correspond to the move of two particles from the orbitals si , S2 
to other orbitals /ix , A*2 ■ At the same time, the rest of particles ( n — 2 particles) can occupy different m — 4 
orbitals. Therefore, the same matrix element stands for other basis states \i') and with the same move of two 
particles, thus, Hyy — = V Sl ^ lS2 ^ 2 . As a result, among matrix elements of the Hamiltonian matrix there 
are equal matrix elements, although they are chosen randomly from the ensemble of two-body random matrices. This 
non-trivial fact indicates that, in spite of completely random character of the two-body interaction, the (non-zero) 
matrix elements of the many-body Hamiltonian are not completely independent variables! 
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FIG. 3. Distribution of normalized off-diagonal matrix elements of the matrix Hij for the parameters of Figs. 1-2. Each 
of matrix element Hij has been normalized to its variance computed from the average over N g — 16 matrices with different 
realizations of the (random) two-body matrix elements. The smooth curve is the best fit to the Gauss which is expected for 
uncorrelated matrix elements. 
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The role of the above underlying correlations which are due to a two-body character of the interaction, is an 
interesting and important problem (see JD|), we will discuss some results in Section 2.6. Here, we would like to show 
that these correlations can be easily detected by the study of the distribution of matrix elements . For this, let 
us take the ensemble of the Hamiltonians H with different two-body random matrix elements V SlS2SaSi keeping all 
other parameters. Then, for any fixed values i and j , we can find numerically the distribution of non-zero matrix 
elements Hij . Finally, we normalize each of these distributions to their variances and make the summation. The 
resulting normalized distribution is shown in Fig. 3. The envelope of this distribution looks like the Gaussian, however, 
the deviations are non- statistical ones which can be easily seen by the x 2 -test (for some bins of the histogram, the 
difference is more than 100 standard deviations). 



B. Density of states and spectrum statistics 

As was pointed out, the density of states p(E) for the two-body random model was found to have the gaussian 
form @|J. Ri gorous proof is given for the limit case of a very large number of particles and orbitals. However, 
even for relatively small values of m and n the distribution is very close to the Gaussian, see Fig. 4. It is known 
that the density of realistic physical systems such as complex atoms and heavy nuclei can be approximated as 
p(E) ~ exp(A \/E — Eq) where Eq is the ground energy. Therefore, the simplest model of the two-body interaction 
does not give exact correspondence to a real density, however, it reproduces very fast increase of the density with 
the energy. One should note that for full random matrices of the standard RMT, the density has the semicircle form 
which is very far from the reality. It is clear that in order to compare statistical properties of our Hamiltonian Hij 
with those of complex quantum systems, one should use the left part of the energy spectrum since the decrease of the 
density in the right part is due to artificial cut-off of the single-particle spectrum (finite values of m ) . 




FIG. 4. Density of states of the TBRI-model with the same parameters as in Figs. 1-3. The average over N g = 20 is taken. 
The smooth curve is the best fit to the Gaussian, with a — 5.72 and E = 25.1 . 



The interesting question is about the type of fluctuations in the energy spectrum of the TBRI-model in comparison 
with those predicted by the RMT. A particular interest is the distribution of normalized spacings S between nearest 
energy levels. For standard random matrices the RMT reveals specific kind of the distribution known as the Wigner- 
Dyson distribution. In particular, for symmetric real random matrices (the so-called, Gaussian Orthogonal Ensemble 
of random matrices), the level spacing distribution P(S) with a high accuracy can be described by the following form, 

P(S) = i7rSexp(~) (12) 

] it was shown that for the limit case of a very strong interaction, when one can neglect the 
d (or, the same, without the leading diagonal in the Hamiltonian Hij ), the form of P(s) 



where (S) = 1 . In ggy 
influence of the mean fie 



8 



turns out to be quite close to the expression (|l2|). On the other hand, when studying the distribution of spacings 
between the levels Ei and with k > 1 , the difference between the result of the RMT and the TBRI-ensemble is 
noticeable and increases with k . This means that the level spacing distribution P(s) is quite insensitive quantity of 
spectral correlations and does not "feel " the difference of TBRI-matrices from the full random matrices. 

In physical applications, the interaction V is typically of the same order as the "unperturbed " Hamiltonian Hq 
since in the mean field approximation the term Hq absorbs regular part of the interaction and V is the (chaotic) 
part of the interaction which can not be included in Hq . However, there are many cases when the interaction is 
weak compared to Hq , therefore, the important question is how spectral fluctuations, in particular, the level spacing 
distribution, depend on the interaction and on total energy Ei . The origin of the Wigner-Dyson distribution ( |l2] ) 
is related to the onset of chaos in the exact eigenstates (see, for example, 0]). In standard random matrices the 
WD-distribution occurs for any energy since all eigenstates arc completely random (their components are distributed 
according to the gaussian distribution for large size of the matrices) . 

Experimental data for complex atoms p3[ , p4[ ] and heavy nuclei |25| (see also references in ||) agree with the 
Wigner -Dyson statistics. The WD-distribution has been also observed in numerical calculations for the Ce atom 
(Hfgll and the nuclear shell-model H-|f]. 
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FIG. 5. Level spacing distribution for the parameters of Figs. 1-3 with N g — 10 . All but 10 levels from both edges of the 
energy spectrum are taken into account, with proper rescaling to local mean level spacings. Solid curve it the Wigner-Dyson 
distribution and the dotted curve is the Poisson, the latter occurs for uncorrelated energy sequences. 

In the TBRI-model, the randomness of the components of eigenstates is different for different eigenstates. This 
fact influences the spectrum statistics if one average over a large part on energy spectrum. In Fig. 5 the level spacing 
distribution for the TBRI-model is shown for the parameters which may be compared with the Ce atom thoroughly 
studied in |l{j-|2jj|, n = 4, m = 11, d = 0.5, Vq = 0.12 . The average over N g = 10 different matrices H+j has 
been done in order to have representative statistics. All but 10 levels from the spectrum edges have been taken into 
account. By searching the structure of eigenstates, one can see that all eigenstates seem to be quite random (see 
examples below). However, one can detect clear non-statistical (regular) deviation from the WD-distribution which 
should be treated as quite strong, having in mind insensitive character of P(s) . This result may be regarded as an 
indication of not strong enough interaction between particles. Detailed study of the onset of the WD-distribution in 
the TBRI model has been recently performed in Ref. p9|] . 

C. Structure of exact eigenstates 

Our main interest in this Section is in the structure of exact eigenstates. The choice of the unperturbed basis, 
reordered in increasing energy, allows us to understand what happens with an increase of the interaction. In the 
standard perturbation theory the natural parameter which controls the intensity of the perturbation is the ratio A of 
the interaction to the mean level spacing D between the unperturbed energy levels. Since the value of D is defined 
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by the total density of states, D — (E) , one can expect that our control parameter is A = Vq/D . However, in 
Ipfjf it was shown that due to a two-body character of interaction, the correct parameter is different. Indeed, in the 
first order of the perturbation theory, the interaction couples not all unperturbed states but those basis states which 
correspond to the shift of not more than two particles. This means that the density of states which are coupled by 
the two-body interaction is much less than the total density p (E) . Therefore, the correct mean level spacing df that 
has to be compared with the interaction Vq , is much larger than D . 

If the interaction is very weak, Vq <SC df , the standard perturbation theory can be applied. In this case any of the 
eigenstates in the unperturbed basis has the form of the delta-function (originated from the unperturbed state |no)) 
plus small admixture of other components with amplitudes decreasing as |C„| ~ 1/ \n — uq\ , therefore, the number 
of principal components is small, N pc ~ 1 . In this case one can speak about perturbative localization of eigenstates 
in the unperturbed basis. This situation is quite typical for eigenstates corresponding to low energies, in the energy 
region where the density of states is small. 

With an increase of perturbation (or when passing to higher eigenstates for the fixed Vq ), the number N pc of 
principal components with essentially large amplitudes C n increases and can be very large, N pc 3> 1 , even if Vq is less 
than df . Such a situation occurs when Vq ^ ^JD df , see details in |1| . In such a case the structure of eigenstates 
is "chaotic " , however, there are many "holes " inside such eigenstates in a given basis. Therefore, in spite of large 
number of components, these sparse eigenstates are non-ergodic, which leads to non-gaussian statistics. Namely, the 
fluctuations of the components C n can be extremely large and statistical description is not valid. One should stress 
that in this case the number of principal components can not be estimated as N pc ss T/D , as is typically assumed in 
the literature (here V is an effective "size " of the eigenstates in unperturbed energy representation, see below). 

When the interaction is relatively strong, Vo ~ df , specific transition occurs from non-ergodic to ergodic eigenstates. 
This transition has been discovered in |30|| by considering the flow of the energy in the Fock-space of excited states. 
For very large number of particles this transition is sharp and may be compared with the Anderson transition in 
solid state models (see for example details in [ |30|]3l| ). Therefore, the condition Vq > df can be considered as the 
transition to chaos inside compound eigenstates, thus, allowing to describe the model in a statistical way, see below. 
An example of such chaotic eigenstates is given in Fig. 6 for the parameters related to the Ce atom. One can see that 
for large excitation energy (n is the number of exact eigenstates reordered in increasing energy ), the eigenstates 
look more extended (delocalized) in the unperturbed basis. It is interesting to note that they look very similar to the 
eigenstates of the Ce atom, obtained in the direct quantum computation based on the Hartree-Fock method [|l9]-^l| . 
For the first time, chaotic structure of eigenstates of the Ce atom has been revealed in p^|. 




FIG. 6. Examples of the exact eigenstates of a matrix Hij for the same parameters as in the previous figures. Components 
Cj of 8 eigenstates for exact eigenstates in the low part of the spectrum are shown. 



Let us now discuss how to quantitatively characterize chaotic eigenstates. First, we introduce the matrix 
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l \ 2 - - 2 



Cj(EW) (13) 



constructed from the exact eigenstates |n) corresponding to the energy £W. In what follows, we use the notations 
which refer low indices to the basis states, and upper indices to the exact (compound) eigenstates. Thus, the structure 

(n) 

of eigenstates is given by the dependence u>j on j for fixed values of n . On the other hand, if we fix the index j 

and explore the dependence Wj on n , one can understand how the unperturbed state \j) is coupled to other basis 
states due to the interaction. The latter quantity is very important since it gives the information about the spread 
of the energy, initially concentrated in a specific basis state \j) , when switching on the interaction. The envelope of 

this function wj™' in the energy representation is known as strength function or local spectral density of states (LDOS) 
and will be discussed in details in next Section. 

From Fig. 6 one can conclude that when the number of principal components is large, N pc 3> 1, such eigenstates 

(n) 

may be treated as random superposition of components C- , although they do not occupy the whole unperturbed 

basis. The gaussian character of the fluctuations of Cj depending on the indices j or n , has been revealed in [l9|-|2l|l 
for the Ce atom, thus allowing to treat the exact eigenstates of a dynamical systems, as chaotic eigenstates. 

The size of the basis which they occupy can be associated with the "size " of eigenstates, or, with the localization 
length. As is known, the notion of the localization length is very important in solid state applications, when studying 
the eigenstates of disordered models in infinite basis in the position representation. In such applications the localization 
length loo is defined via exponential decrease of the square of the amplitude of eigenstates . For the finite basis, the 
definition of the localization length can be generalized in a way described in H] . Following to |7]]|] we define here two 
localization lengths, namely, the "entropy localization length" lh, and the localization length li pr associated with the 
so-called participation ratio. The first one is defined by the expression 

l h = Nexp((H)-H ), (14) 

where (Ti) is the mean " entropy" of eigenstates, 

M N 

< w > = -^EE w 5" )ln (^ n) ) ( 15 ) 

n=l j=l 

and Tio is the normalization constant which is equal approximately to 2.08 in the case of pure gaussian fluctuations 
of Cj (see details e.g. in 0). Here M is the number of eigenstates which are taken for the average. This can be the 
number of eigenstates of one matrix Hij taken from a small energy window, or the number of eigenstates for the fixed 
n , computed from different matrices with different realization of disorder in two-body matrix elements. 

The second definition is commonly used in solid state applications (see e.g. H]). Assuming the gaussian character 
of fluctuations of the components of eigenstates, the localization length li pr is defined by 

At N 

V = |; ^EE(-r) (16) 

n=l j=l 

In the above expressions (|lj, |l6|) the factors 2.08 and 3 are, in fact, normalizing coefficients which provide, in the 
limit case of completely extended and gaussian eigenstates in the finite basis of the size N , the "maximal" value of 
the localization length lh = k pr = N. In other extreme case of a strong (exponential) localization in the unperturbed 
basis, the above two localization lengths are proportional to that found from the tails of eigenstates (see 0). The 
dependence of the localization lengths lh and h pr on the basis number and energy of exact eigenstates is given in 
Fig.7a-b. The average over N g — 50 matrices with different realizations of two-body matrix elements has been 
taken, in order to smooth strong fluctuations in the value of localization lengths of individual eigenstates. The above 
definitions of localization length can be taken for the estimate of a number of principal components N pc , and can be 
associated with the degree of "chaoticity " of compound eigenstates. The data of Fig. 7 show quite good correspondence 
with direct computations fl9| , ^o]| performed for the Ce atom, for which the localization length lh was found to be 
about lh ~ 110 — 130. Comparing Fig. 7b with Fig. 4 of the density of states p (E) , one can see the similarity. This 
fact can be understood from the simplest estimate of N pc as a total number of (basis) states defined by the "width" 



r , the latter can be approximated as the mean-square-root of the distribution Wj for the fixed 

r 



N pc ^-=Tp(E) (17) 

One should stress that this expression is valid for ergodic eigenstates ( Vb > df ) and shows the proportionality of 
N pc to the density. As for the width T , it is approximately independent of the excitation energy, see below. 
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FIG. 7. (a) Localization lengths 1^ and U pr are given in dependence on the number n of exact eigenstates |n) , see ([Tij) and 
(p^[). The parameters are taken the same as in Figs. 1-6, with N g — 50 . (b) the same as in (a), but in the energy representation. 



Eq. ([i5|) reflects general relation between the number of principal components N pc and the entropy of eigenstates, 



S 



EF 



In N„. 



(18) 



where Sef stands for any reasonable definition of the entropy. In application to shell models of complex nuclei this 
relation has been studied in great details in @,^8| . Combining Eq. ([18|) with Eq. ([l7|) , one can get 



Sef « In 9 (E) + In T(E) 



(19) 



In contrast to the density of states, the width T is a weakly dependent function of the energy. Therefore, the 
entropy Sef found from exact eigenstates practically coincides with the thermodynamical entropy, the fact which 
was mentioned for the first time in |E6L| 



D. Strength function 

In this Section we discuss the properties of the strength function which is defined as 

W{E {m \ j) = \cl n) \ 2 S(E - E&>) (20) 

71 

Here C J ( ™ ) are the components of eigenstates \n) ("compound" states) of the total Hamiltonian H given in the 
unperturbed basis \j) , and E^ is the energy associated with the state \n). The sum is taken over a number of 
eigenstates \n) chosen from a small energy window centered at the energy . One can see that this function 

W(E, j) is originated from the same matrix which has been introduced in previous Section when discussing the 
structure of exact eigenstates. Indeed, an exact eigenstate is characterized by the dependence on j for the fixed 
value n (associated with the energy E^ n ' ). On the other hand, the strength function is characterized by the same 
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function when index j is fixed and we are interested in the dependence on the energy ij( m ) due to the relation 
between 2?( n ) and n therefore, 



W{E^ m \j) ~Fj n) p(E), F^ n) =wf\ (21) 

Here we have introduced the F— junction Fj which gives the envelope of in dependence on the indices j and n 
(the bar stands for the average inside small windows centered at j and n). In fact, the strength function is the (smooth) 
representation of a (simple) basis state \j) in terms of exact eigenstates. This function is very important since it can 
be measured experimentally. It contains an information about the internal interaction between unperturbed states. 
Namely, it shows how the unperturbed state \j) is coupled to the exact states \n) due to the interaction. An effective 
width of this function (spreading width ) defines the energy range associated with the "life time" of an unperturbed 
state \n) if initially one excites specific basis state. 

In solid state models the role of the unperturbed energy in Eq.(|2^) plays the position j of an electron and the 
function W(E^ m \j) =^ W(E,j) has the meaning of the electron density of states for the fixed position j . For 
this reason this function is known in solid state physics as the local density of states (LDOS). One can see that if 
(apart from fluctuations) Wj is independent of the position (or, in our application, the energy of compound state), 
it reduces to the total density of states, W (E^ m \ j) p (E) . Also, if the total density of states is constant, p = a -1 
, the dependence of the LDOS on j is of the form W(E, n) = W(E — aj) . This means that the form of the LDOS is 
the same for any j . In our case of strong dependence of the density of states on the energy, the form of the LDOS is 
quite complicated, see next Section. Normalized to the mean energy level spacing, the strength function W(E^ m \ j) 
determines an effective number N pc of principal components of compound states \n) which are present in the basis 
state \ j). 

Similar to the analysis of the structure of the eigenstates in dependence on the interaction, one can understand 
that for a very weak interaction Vq <C df the LDOS is a delta- like function with a very small admixture of other 
components which can be found by the standard perturbation theory. With an increase of the interaction, the number 
of principal components increases and can be very large. However, if the interaction is not strong enough |i~5| ] , 

1»^»±M, (22) 

d f ir 2 y d f 

the LDOS is sparsed, with extremely large fluctuations of components, see details in Ref. pi] ]. In order to have ergodic 
LDOS, one needs to have the perturbation large enough, Vq 3> df (for a large number of particles this transition is 
sharp and, in fact, one needs the weaker condition, Vq > df , see details in [151). 




FIG. 8. Few examples of the LDOS for the TBRI-model for different basis states j = N/8, JV/4, 3JV/4, N/2 . The average 
is done over a number of eigenstates with close energies, for n = 4 particles, m = 11 orbitals, Vq = 0.12, do = 0.5 and with 
additional average over iV g = 50 matrices Hij . Smooth curves are the best fit to the Gaussian at the center of the energy 
spectrum, j = N/2 . 
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For Vo 3* df the LDOS turns out to be ergodic and thransition to chaos occurs J30| , p3| , p 1|J34|J29| , |35 37|, therefore, 
statistical description of the model is valid. Few examples of the LDOS for the TBRI-model are given in Fig. 8 for the 
parameters of the Ce atom. One can see strong dependence of the shape of the LDOS on the position of basis state 
|j) in the energy spectrum. 

Before we start with the discussion of analytical results, it is important to point out the correspondence between 

the shape of the LDOS, and that of exact eigenstates which have been discussed in previous Section. Specifically, 

(n) 

from the analysis of the structure of the eigenstate matrix ur- one can expect the similarity between the LDOS 
and shape of the eigenstates. For the first time such a similarity has been observed when studying band random 
matrix ensembles |^,^9|. Moreover, the detailed study of some dynamical models [ p"7p0| ] have revealed that even 
when the shapes of the LDOS and EFs seem to be completely different, after a proper rescaling which involves both 
unperturbed and perturbed energy spectrum, both shapes are very similar. In order to compare characteristics of the 
LDOS and EFs, in Fig. 9 the dependence of the entropy localization length on the energy is shown for both the LDOS 
and EFs. Apart from strong fluctuations (which are due to a chaotic nature of the components ), in general, the 
dependencies Ih (E) look very similar. 

The problem of the correspondence between shapes of the LDOS and EFs is still open, however, from the studies 
made up to now (see also J17,[Io||), one can conclude that if the width of the perturbed spectrum is of the same order 
as the unperturbed one (or, the same, the perturbation is not very larger), one can expect that both shapes are very 
close to each other. The importance of this problem of similarity for the shapes of the LDOS and EFs will be clear 
in next Section when we discuss the relation between the shape of the EF and generic properties of the occupation 
number distribution. 




For the first time the form of the LDOS in random matrix theory has been discussed in Ref. JoJ where band random 
matrix ensemble has been considered. The matrices were assumed to consist of the diagonal part Hq in the form of 
reordered numbers Hjj = j D (thus, the unperturbed level density is constant, p = D^ 1 ), and the perturbation Vq 
with the random independent off-diagonal matrix elements for \i — j\ < b. Outside the band of size b 3> 1 matrix 
elements are zeros. The distribution of the off-diagonal matrix elements is characterized by the zero mean, (Hy) = 
and the variance V 2 = {Hfj) ■ For this model, the relevant parameter was found to be 

q=^±. (23) 

Wigner analytically proved || that for relatively strong perturbation, Vq ^ D in the limit q <C 1 (when the influence 
of the main diagonal Hjj is strong) the form of the LDOS is the Lorentzian, 

Wbw (-E) = j- - , E = E-Dn (24) 

Z7T _|_ 1 BW 

which is nowadays known as the Breit-Wigner (BW) form. Here Tbw is the spreading width which is half- width of 
the distribution (pj), 
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BW 



2ir Po V* (25) 



and the energy E = E — D n refers to the center of the distribution. In other limit q 1 the influence of the 
unperturbed part Hq can be neglected and the shape of the LDOS tends to the shape of the total density of states 
of band random matrices without the leading diagonal Hjj, which is the semicircle. 

Recently, the shape of the LDOS has been studied rigorously for a more general distribution of the off-diagonal 
elements v nm in the Wigner band random matrix model (WBRM). Namely, the variance of random matrix elements 
is taken to depend on the distance r = \n — m\ from the principal diagonal according to an envelope function f(r) 



which, for r — > oo, decreases sufficiently fast (see details in 41 ). In this case the effective band size b is defined by 
the second moment of the function f(r) . Another important generalization of Ref. |4lJ is related to the sparsity 
inside the band, which can be defined as a relative number of zero elements in each line of the matrix. As was shown 
above, such a sparsity, which is due to the two-body nature of the interaction between the particles, is important for 
statistical properties of compound states. 

The BW-form of the LDOS also occurs in a more general model describing the interaction of an unperturbed state 



with a large set of complex states, see details in the book |42|. It also appears for dynamical systems with complex 



behavior, such as the Ce atom ]l9|-|2l[|, the sd shell model p6|-28|, as well as in random models in application to solid 
state physics [f|||J]. 

For a long time it was believed that the LDOS for complex physical systems has the universal energy dependence 
described by the BW-distribution (pi|). However, when studying the structure of the LDOS and EFs of the Ce atom, 
it was observed some deviation from the BW-shape for large distance from the center. By applying the WBRM, in 
JlH it was analytically shown that the shape of the LDOS is highly non-universal for energies larger than the effective 
size of the band in the energy representation, \E\ ^> Db (in what follows by E we mean the distance from the center 
E c of the distribution, E =>• E ), see also numerical data in |^|. Namely, outside this range, the tails of the LDOS 
are highly non-trivial, decaying very fast (even faster than the exponent) when \E\ — > oo . Therefore, the range of 
parameters for which the form of LDOS has the BW-form in the WBRM-model, is given by the condition [ElJ which 
can be written as 

Po 1 < r BW < b Po 1 (26) 

Here, the left-hand-side of the inequality is related to non-perturbative character of the coupling since the perturbation 
should couple many unperturbed states. 

The condition (Kq) turns out to be of generic and can be applied to real physical systems. Moreover, it can be used 
to find the effective band- width b e ff of interaction when other definitions of b are obscure. In the TBRI-model with 
many interacting particles the band- width is very large (practically it is infinite), and does not play any role in many 
solid state applications. On the contrary, in application to complex atoms and nuclei, the number of particles above 
the Fermi level is relatively small (four particles in the Ce atom and 12 particles in the sd shell models), and this 
effect can be important. It should be pointed out that the BW-dependence has infinite second moment. At the same 
time, in any physical application the second moment is always finite and is defined by the sum of the square of the 
off-diagonal elements. This fact results in non-Lorentzian tails of the LDOS. 

The important question is about the shape of the LDOS for a strong interaction, when the condition ( |26] ) violates, 
r bw ~ bpQ 1 . Since in the mean field approximation the "regular " part of the interaction is included in the mean 
field, very often the interaction V is compared with the unperturbed part Hq ■ Therefore, this situation when the 
shape of the LDOS is very different from the BW-form, is quite physical. In order to understand what happens in 
this case, it is useful to study how the transition from the BW-form to the semicircle occurs in the WBRM-model. 
As was pointed out, the semicircle form itself is unphysical and appears when neglecting the unperturbed part Hq . 
However, the form of the LDOS in the transition region seems to be of quite generic. 

Numerical data in Refs. pl[ , }45]] for the WBRM-model have shown that for the case when Tbw ~ bpg 1 the form 
of the LDOS can be approximately described by the Gaussian. The same fact has been observed and discussed in 
Refs. [ p6|j28| when studying eigenstates and LDOS for shell model of nuclei. Since in the latter application the band 
size of the interaction is not well defined, it is better to introduce one more parameter, in addition to the half-width 
Tbw ■ As was shown in Ref. p5[ ], in general, the form of the LDOS can be effectively described by two independent 
parameters, Tbw and AE where the latter is defined via the variance of the LDOS, 

(AEf = 4a w = ij(E- E c f W BW (E) dE (27) 

In the symmetric case (in the TBRI-model, at the center of the spectrum) , E c coincides with the unperturbed energy, 
E = E j . 
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Then, if the parameter Tbw is much less than AE , the shape W (E) of the LDOS is the BW-dependence ( |24| ) and 
Tbw has the meaning of the half-width of the distribution W(E) — Wbw (E) . On the other hand, if Tbw ~ AE , 
the form of the LDOS is approximately the Gaussian, 



W{E) 



1 



a w 



exp 



(E-Ec. 
2(a w f 



(28) 



Detailed numerical study |45|] of the form of the LDOS in the region Tbw ~ AE have shown that the LDOS coincides 
with the Gaussian with a very high accuracy. It is important to stress that in this case, the parameter Tbw has 
nothing to do with the half-width Thw of the LDOS. the latter is proportional to the mean-square-root aw of the 
Gaussian, Thw w Cq aw (see also discussion in p6| , p8| ). Note that the Gaussian form typically occurs in "statistical 
spectroscopy" ]To| , ^6| , ^7| when neglecting the mean field term Hq in Eq.(|l|). One should note that sometimes the 
fact that the LDOS can deviate from the BW-form due to the influence of the finiteness of AE , is missing in the 
literature. 

It is important that the center and variance of the LDOS can be explicitly expressed via diagonal and off-diagonal 
matrix elements respectively jll|. Indeed, the center is defined by 



e c = =j2 E(n)F j n) 



E 



01 



(n\ H \m) 



\ m \ J, 



= H n - e j 



(29) 



where the relation (i \ H \j) = Sij (i \ H \j) is used for exact eigenstates. Correspondingly, the variance can be obtained 
from the matrix elements of H 2 , 



{°w)j = E ( E(n) - 2 ^ « E ( Ein) - e,) 2 = E * 



2 

jp 



(30) 



For the TBRI-model the sum of the off-diagonal elements for any fixed value of p can be evaluated exactly [fL3[ 



(&EYi 



J2 H j P = v o ( n - + v o K 2 = jV£n(n - 1) (m - n) (3 + m - n) 



(31) 



where the expressions (Q) have been used for K\ and Ki , and Vq is the variance of the off-diagonal matrix elements. 
One can see that with an increase of the interaction, the half- width of the LDOS changes from the quadratic dependence 
Thw ~ Vq to the linear one, T^w ~ Vo ■ K is interesting to note that the variance of the LDOS for Fermi-particles 
does not depend on the index j which stands for a specific basis state, therefore, (AE) 2 (AE) 2 



E. Analytical solution for the LDOS 

Very recently, the form of the LDOS for the TBRI-model has been analytically found in (l6) for any strength of 
perturbation. In this Section we discuss the approach of Jl6| and the obtained results. We would like to stress that 
our aim is to find the LDOS in terms of matrix elements of the total Hamiltonian Hij , without its diagonalization. 
To start with, let us rewrite the general expression for the LDOS in the form 

W k (E)=F(E k ,E)p(E) (32) 

where E is the total energy of the system (energy of an exact eigenstate). As was pointed out, the F— function gives 
the shapes of both exact eigenstates and strength functions depending on what is fixed, the total energy E = E^' 
or the unperturbed one, E k . The method used in [|l6| is an extension of the approach developed in [ 
takes into account specific structure of the Hamiltonian TBRI-matrix. Specifically, first, we fix some basis component 
|fc) and diagonalize the Hamiltonian matrix without this component. Then the problem is reduced to the interaction 
of this component with exact eigenstates \i) which are statistically described by the matrix components V k i- In the 
spirit of the approach of Ref. ]42| ] let us introduce a small energy window A which will be used for an average over the 
total energy inside this interval. As a result, the set of equations for the LDOS can be written in the following form, 



1G 



H ' <£ » - I jE^WTJrm (33) 



where 



T k (E)~2v\V kl \ 2 P {E) (34) 



is some function which can be associated with the half- width of the distribution W k (E) , for the case when the energy 
dependence is weak, T k (E) ~ const . The energy shift S k for the basis state \k) , 

^ \V kl \ 2 (E-E^) 

Sk = Z>iE-EMy + % (35) 

is due to the asymmetry of the perturbation, if the energy E k is not at the center of the spectrum. This shift is just the 
modified second order correction to the unperturbed energy level. For the calculation of the shape of the eigenvector 
\i > one should substitute the exact energy E = 25 W — Ej + S^. Then, if the interaction is not very strong, in the 
evaluation of the above equations the difference — 5 k can be neglected. 

One should stress that the summation in the above equations is performed over exact states. Since the exact 
eigenstates are unknown, one should express everything in terms of the basis states only. To do this, we express exact 
eigenstates \i > through the basis components, 

|V«l a = J2 \ C p } \ 2 \H kp \ 2 +Y. C( : ] * C P )H ^H q u (36) 

P Pt^Q 



In previous Sections it was argued that for large number of principal components, N pc 3> 1 and sufficiently strong 

interaction Vq > df , the components Cp} q can be treated as random variables, therefore, the second term in ^ 
vanishes after averaging. Substitution of Eq. (^6|) into Eqs. (Q [35]) gives 



T k (E) = 27^1^1^(25) = J2 \ H *p\ 2 Tp(E) the) ( 3? ) 

P ^k p^k (E P + dp — E) 2 + — 



5S E - E " tfi(E-E,-S p )' + ^- 

where the integral is taken as the principal value. Last equality is valid in the approximation of slow variation of 
T p (E) and 5 p . The equations for T k (E) and S k allow to calculate the strength function ( (33| ) from the unperturbed 
energy spectrum and matrix elements of the total Hamiltonian H. 

Now we have the set of equations d3^ , |3^ , |38| ) which, in principal, give the solution for the LDOS W k {E) and can 
be solved numerically. However, for relatively large number of particles (practically, for n > 4), one can find an 
approximate analytical solution of the problem |l6[ ]. By analyzing these equations, in Jl6[ it was proved that the 
condition of the self-consistent solution for the LDOS (or the same, when the shape of the LDOS exists as a smooth 
function of the energy), is just the condition for the onset of chaos in the TBRI-model, Vq > df , discussed in previous 
Sections. 

More specifically, for a very strong interaction, L >> df the number Nf of effectively large terms in the sums is 
large, Nf ~ T/df, fluctuations of L are small, 6T ~ Y / ^/Nf and Eq. (|37| ) can be written in the form, 



T k {E)^2n\H k p\ 2 pf{E) (39) 

Here E = E — 5 and the energy shift 5 = < 5 p > can be neglected in the case of L < < a with a standing for an effective 
band- width a of the Hamiltonian matrix H pq (see Eq.(|4l|)). In order to perform the summation over p, it was assumed 
that r(25) and pf (25) change slowly within the energy interval of the size T. As a result, in order to have large number 
of final states Nf ~ 2ir H 2 p /d 2 and statistical equilibrium (small fluctuations of L), one needs H kp >> df. In this case 
chaotic components of exact eigenfunctions in the unperturbed many-particle basis ergodically fill the whole energy 
shell of the width L, with Gaussian fluctuations of the coefficients Cjp and the variance given by the F— function (|2l| ) 



(see also MO* 
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With a decrease of the ratio Hk p /df , the fluctuations of T increase and for Hk p < df the smooth self-consistent 
solution of Eqs.( [37|) d oes not exist. Indeed, in this case the term T p in the denominator of Eg. ( ^7|j3^ ) can be neglected 
and the sum in (B7j) is dominated by one term with the minimal energy E — E p ~ df. Therefore, for a typical 
basis state \k > formally one gets ~ T p (Hk p /df) 2 << T p . This contradicts to the assumption of the equilibrium 
according to which all components are of the same order, ~ T p . 

One should stress again that the absence of a smooth solution for the shape of the eigenstates and the strength 
function does not mean that the number of principal components N pc in exact eigenstates is small. It can be large, 
however, the distribution of the components is not ergodic, there are many "holes" inside exact eigenstates which 

occupy the energy shell of the width 2ir\Hk p \ 2 Pf(E) (see pp| , p^| ). In such a situation, the fluctuations of are 
very large and non-Gaussian. 

It is important to note that the ensemble average (over many matrices U^ p ) in this problem is not equivalent to 
the energy average (inside specific Hamiltonian matrix). Indeed, the average over the single-particle spectrum leads 
to the variation of energy denominators in ( |37|) and can fill the "holes " in the F— function. 

From general equations for the shape of the LDOS one can make an unexpected conclusion that the spreading width 
T(E) can be a strong function of excitation energy E due to the variation of the density pf (E) = dj 1 of final states 
in Eq. (|39|). In Ref. |l6| it was shown that for the excited states, well above the ground state, the energy dependence 
of pf (E) and T(E) can be quite close to the Gaussian. This result is based on the estimate of the two-body density 
Pf (E) , 

p f (E)=pf(E)+pf(E) (40) 

(2) 

where the density pf is determined by the energy difference oj pk ! between the states \p > and \k > which differ by 
the position of two particles, and by lu^ between those states which differ by the position of one particle. Detailed 
analysis have shown that both pf (E) and pf (E) for large number of particles are described by the Gaussian, 



(1,2), _ Kl,2 

P f ( E ) 7JF= exP 

cri i 2v27r 



V 



E-E k - 
2^ 



(41) 



Normalization parameter ifx,2 stands for the number of one or two-particle transitions, see Eq.(^). Here the average 
frequency of one and two-particle transitions reads as 

iiW»m/(m-n)(e-£ii/n) (42) 

and 

^=2(^-^)«2m/(m-n)(e-£ fc /n) (43) 

where ejt = E^jn is the mean single-particle energy in the basis state \ k > containing n particles, e is the single-particle 
energy averaged over all m orbitals, and the mean energy of the empty orbitals can be found from the relation 
me = efcYi + e^(m — n) . 

,(1.2), 



The variance of pf '(E) for one and two-particle transitions is 



a\ = a 2 p + a\ + 2(n - l)V 2 w a\ + 2(n - 1)V 2 (44) 

and 

o\ = 2a 2 p + 2a 2 + (An - Q)V* « 2a 2 + (An - 6)V 2 (45) 

where a 2 is the variance of single-particle spectrum, and Vq is the variance of non-diagonal matrix elements of the 
residual interaction. Note that in the case of n < < m for low- lying states the variance of the occupied orbital energies 
a 2 is small and the variance of empty orbital energies is a 2 ~ a 2 . 
Thus, the width Y(E) is given by the following expression, 

T(E) = 2vr \(n - l)V 2 pf (E) + V 2 pf (E)] (46) 
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The factor n — 1 appears since for single-particle transitions the summation in H kp = 5Z„ V av -*^v is performed over 
occupied orbitals. When the ratio i^ 2 /((n — l)-^i) = (m — n — l)/4 is larger than 1, the two-particle transitions 
dominate and one can neglect the differences in uJ and a for two-particle and one-particle transitions. In this case the 
spreading width is described by the simple Gaussian form, 

T k (E) * 27 r( AE)l -^=e Xp { - &-E k -u*)> 1 

where E = E — 5. Here (A.E)? is the variance defined by Eq.(|3l]) and Wfe and a k are close to that for the two-particle 
transitions. The maximum of pf (E) and T(E) is shifted by the value \tJk\ towards the center of the spectrum, 
compared to the maximum of the Breit-Wigner function. This leads to some distortion of the strength function 
Eq.(|33|) and the shape of the eigenstates, which is especially large at the bottom of the spectrum. 

0.0025 I 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

W k (E) - 




k 



FIG. 10. Shape of the LDOS ( p3| ) in the basis representation. Broken line is the result of numerical diagonalization of the 
TBRI-matrix. To reduce fluctuations, the average over 50 different matrices Hik and over a number of nearby components has 
been made. The computation has been done for n = 6 particles, m = 13 orbitals, therefore, the total size of the matrix is 
N = 1716. The interaction strength is Vq ~ 0.1 and do = 1-0 . Dashed and smooth full curves obtained by computation of 
Eq.@ with T k {E) given by Eqs.(H|i|) and by Eq.@ correspondingly. 
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FIG. 11. Dependence FfJE) is shown for the parameters of Fig. 10. Full curve is the expression (^), the dashed curve is 
the computation from Eqs.(p7], pq ). 

The above consideration shows that if the interaction is small, V <C <Jk , the strength function has the Breit-Wigner 
shape with a broad gaussian envelope described by the dependence Tk(E) in the numerator of Eq. (|33|). In fact, such 
a dependence results in the correct (finite) variance of the strength function. When the interaction V increases one 
needs to take into account one more contribution in Eq.p^) (it was neglected in Eq.(|39|)). It increases the width of 
W P (E) and leads to the estimate a\ ~ <j\ + V*. With further increase of interaction, where the shape of W P (E) is 
close to the Gaussian, one gets o\ ~ erf + (AE)^. 

Direct numerical study of the model (Q) with n — 6 Fermi-particles and m = 13 orbitals confirmed that the 
above analytical expressions give quite good description of the shape of the strength function Wk (E) as well as the 
energy dependence Tk(E) , see Fig. 10 and Fig. 11. The size of the Hamiltonian matrix is N = = 1716 and the 
unperturbed state io = 440 was taken. 

For the case of quite strong interaction, when r ~ a, the gaussian variation of T(E) in the numerator of Eq.(|3^) 
becomes as important as the variation of the Breit-Wigner energy denominator (E — Ek) 2 + (r/2) 2 . In this case the 
transition from the Breit-Wigner type to Gaussian shape of the LDOS (strength function) takes place. However, still 
one can use Eqs.(|37[ |3|) and ( p3| ) in order to calculate numerically T(E), Pk(E) and F(E, E k ), taking T from Eq.(f47|) 
with a\ ~ a\ + (AEyj^as the zero approximation in the right-hand side of Eqs.(|37|, |38| ). 

F. Non-statistical properties of the TBRI-model 

In previous Sections we have considered statistical properties of the TBRI-model based on chaotic structure of the 
eigenstates and LDOS. However, one should be very careful with this approach since due to underlying correlations 
in matrix elements of H (see Section 2.1.3), the approach is not always valid, even if all two-body random matrix 
elements are completely random and independent variables. This fact is due to a two-body nature of interaction and 
should be taken into account in some cases. Below we show an example when statistical description is incorrect (see 
details in @). 

Let us consider a single-particle operator 

a.fi a. .13 

where aj, and ap are the creation and annihilation operators and we have introduced the density matrix operator 
Pap = a„fl/3 which transfers a particle from the orbital f3 to the orbital a. The matrix element of M between compound 
states can be expressed through the projection of the density matrix into the basis states, 

(m| M \n 2 ) = M af3 (m| PaP \n 2 ) = ^ M aP P { ^' n2) (49) 

a/3 a/3 

where 

p^^^)^^^ (5Q) 

ij 

is determined by the exact eigenstates only. In what follows, we are interested in statistical properties of this many- 
body operator p(a, (3) for the fixed orbitals a and /? which, on the other hand, determines statistical properties of 
the single-particle operator M , see Eq.(f48|). 
One can see that this operator has zero mean, 



p(a,j3) = (m| Pa/3 |n 2 ) = (51) 



if compound eigenstates are truly random. 

In general case, the variance of p(a,[3) which is of our main interest, has the form, 



p 2 (a,/3) = {ni\p a p \n 2 ) (n 2 \p0 a K 
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C^C^C^C^ (i\ Paf3 \k) (l\ P(ia \j) = s^ n2) + si n ^ 

k,l 



Here we separated the diagonal, 



S 



(n 1 n 2 ) 



(ni) 



ik 



c. 



\(i\p a p\k}\ , 



and non-diagonal, 



(52) 



(53) 



(54) 



contributions to the sum (|52j) and assumed that eigenstates are real vectors (note that our matrix is symmetric) . 
Typical shape of the density matrix is shown in Fig. 12. This shape can be compared with the statistical approach 
developed in Ref. J4(|[l7| for a very large interaction, in the case when the role of the unperturbed part H is neglected 
(therefore, the influence of the leading diagonal Hjj is small). 




FIG. 12. Form of the density matrix (pOI) for the TBRI-model with the parameters of Fig. 1-5. The quantity 

1/2 



Pn,m 



f (n 1 n 2 ) 
\Poc0 



is computed with the average inside the blocks of size 10 x 10 (in the same way as in Fig. 2). 



As we have discussed in previous Sections, for sufficiently large interaction compound eigenstates \n) of the TBRI- 

(n) 

model may be considered as pseudo-random functions due to a very large number of components C\ . Therefore, 
it is natural to expect that the non-diagonal part S^" 1 '" 2 ' is zero and the variance is essentially determined by the 
diagonal term g^™ 1 '™ 2 ) which can be described statistically (for this statistical approach see Refs. @|19||0|1). This 
assumption has been used in the previous calculations of matrix elements between compound states in ] 49[ 19|, |50| , |5l|] . 
However, recently it was shown |13|] that in many-body system these two terms are of the same order, S c ~ Sd, 
even for completely random two-body interaction. In order to show this very unexpected effect, in Rcf. [ fl3[ direct 
computations of the terms S c and Sd have been performed for the TBRI-model, see Figs. 13-14. One can see that 
the data reveal a systematic difference between the diagonal approximation and exact expression ( |52| ) . In particular, 
Fig. 14 shows that non-diagonal term S c is of the same order as Sd which clearly indicates the presence of correlations. 

Below, we show how these correlations emerge in the non-diagonal term S c (for more details see Ref. |0|). First, 
note that for a given i the sum over k in Eq. ( |53| ) for Sd contains only one term, for which |A;) = a^a a \i) = \i'), 
determined by transferring one particle from the orbital a to the orbital (3 in the state \i) (hereafter we use the 
notation i' to mark such states). Accordingly, the index i runs over those states in which a is occupied and (3 is 
vacant. For such i and i' the matrix element (i\ p a p — 1, otherwise, it is zero. Therefore, in fact, the sum in ( [53] ) 
is a single sum, with a number of items less than AT, 
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-i(n 1 n 2 ) 



(55) 



where the sum runs over the specified z. Analogously, Eq. (54) can be written as the double sum over i and j 
specified as above, 



j(nins) _ p("i)p(ni)gr(i2)p(n2) 
c / ; 2 j i j' 5 



(56) 



where j' is a function of j, = dga a \j). Note that the energies of the basis states and their primed partners are 
connected as Ei> — Ei = ep — e a = Ey — Ej. 



0.001 




3000 3000 

FIG. 13. Mean square matrix element ( |52| ) calculated in the TBRI-model for n = 7 particles and m = 14 orbitals, a = 7, 
f3 = 8, as a function of the eigenstate 712 for ni = 575 (total size of the matrix is N = 3432). Dots correspond to the sum 
5d + S c while the solid line represents the diagonal contribution Sd only [see (p3J)] . 




1000 2000 3000 

FIG. 14. Ratio R = S c /Sd of the correlation contribution to the diagonal contribution for the same parameters as in Fig. 13. 
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One can expect that maximal values of the sum ( |55| ) and ( |56| ) correspond to the terms for which C— components 

are principal components of the eigenstates. This means that mean square of the matrix element |(rai|p Q /3 ^2) | 2 
is maximal when the operator p a p couples the principal components of the state \ni) with those of \n 2 ), i.e. for 
jr;( n i) — E( n 2) UJa ^ = £q _ p ar f rom the maximum (\E^ ni ' > — £;(" 2 ) — u) a p\ > T) a principal component of one 
state, say, n\, is coupled to a small component k of the other state n 2 (|-E/c — £ ,< - Tl2 - ) | > T). The latter case is simpler 
to consider analytically, since the admixture of a small component in the eigenstate can be found via the standard 
perturbation theory. This approach reveals the origin of the correlations in the sum S c , Eq. (|5^). For example, if 
cj™ 1 ^ is a small component of the eigenstate ni, then it can be expressed as a perturbation theory admixture to the 
principle components. If the component C^' 



is one of the latter, then there is a term in the sum (pf 

2 



proportional to the principal component squared, 



a 



(«i) 



Based on this consideration, in Ref. 
non-diagonal terms read 

q(n 1 n 2 ) 
1 j ,. 



L3l was found that far from the maximum, 



which is 



> T, the 



A similar calculation of the diagonal sum Sf in2 \ Eq. flU, yields 





2 











(57) 



S. 




("i) 



(1112) 
d 





2 











H 2 



(58) 



Let us estimate the relative magnitudes of Sd and S c 
= "^a^a^a^i), in this case H tj = V„, 



orbitals, \j) = </;, </,,, »;..« 



First, we consider the case when \i) and \j) differ by two 
Since the basis states \i'} and \ f) must differ by the same 



two orbitals, we have Hi<y — V Ulflll , 2tl2 = Hij (note that vi, fi\, v 2 , M2 7^ &,(3, since both states \i) and |j) contain a 
and do not contain /?, whereas and \f) contain (3 and do not contain a). Therefore, the averages over the non-zero 
matrix elements between such pairs of states give HijHi>j' = if?- = Hf,j, = Vq. 



Now, let us consider the case when \i) and | j) differ by one orbital | j) 



= «t 



only. In this case the Hamiltonian 



matrix elements are sums of the n— 1 two-body matrix elements, see Eqs.(^) and (|l0|). As was shown in Section 2.1.3, 
the sums of n — 2 terms in and H^ji coincide and the difference is due to the one term only (orbital a is replaced 
by the orbital (3). Thus, 



HijHyj, = (n - 2)V Z 



{ H i]) 2 - {Hi'j> f 



{n-\)V 2 



(59) 



where we took into account that V K Xf_n,V Kl \ llllUl = V l 8 KKl 5\\ l 5 l _ lt _ ll 5 UUl . 

The contributions of one-particle and two-particle transitions in Eqs. ( |57| ) and (^) representing S c and S4 respec- 
tively, are determined by the numbers of such transitions allowed by the corresponding sums. For the single-prime 
sums in Eq. (|5^) these numbers are proportional to K\ and K 2 , Eq. (Q). In the double-prime sum in Eq. (|57j) 
these numbers are proportional to K\ and K 2, the numbers of the two-body and one-body transitions i — ► j, in 
the situation when one particle and the two orbitals (a and f3) do not participate in the transitions. These num- 
bers can be obtained from Eq. (0) if we replace n by n — 1, and m by m — 2, so that i^i = (ra — l)(m — n — 1), 
K 2 = (n— l)(n — 2)(m — n— l)(m— n — 2)/4. Finally one can obtain that for li^" -2 ) — E^™ 1 ) — > T the contribution 
of the correlation term to the variance of the matrix elements of p a p can be estimated in the ratio as 



R- Sc 

R= s~ d 



(n - 2)K 1 + K 2 
(n-l)K 1 +K 2 



(n — 2)(m — n — l)(m — n + 2) 
n(m — n) (m — n + 3) 



(60) 



This equation shows that for n — 2 we have S c — 0, which is easy to check directly, since in this case Hi'jiHij = . 
For n > 2 the correlation contribution S c is negative in the tails of the LDOS. This means that it indeed suppresses 
the transition amplitudes off-resonance (see Figs. 13-14). For n, m — n ^> 1 the ratio R approaches its limit value —1. 
It is easy to obtain from Eq. (60) that for m — n 3> 1 
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Sd + s c 
s d 



= 1 + R: 



2m 

n(m — n) 



(61) 



Thus, surprisingly, the role of the correlation contribution increases with the number of particles. This result is 
supported by numerical data reported in Ref. Jl3| . In Fig. 14 one can see that the suppression of the matrix elements 
p 2 (a,(3) due to the correlation term at the tails, is quite strong, numerical ratio is R « —0.7 vs. R = —0.55 
obtained from Eq.(|60|); this should be compared with the case n — 4, to = 11, JV = 330 for which numerical value is 
R ~ —0.45, see details in jF3|. The correlation contribution should be even more important in compound nuclei, where 
N ~ 10 5 . This case can be modeled by the parameters n = 10, m = 20; then we have R = —0.66, or, equivalently, 
(Sd + S c )/Sd — 0.34, which means that the correlations suppress the squared element p 2 between compound states 
by a factor of 3 (far from its maximum). 

It is worth emphasizing that the existence of correlations due to the perturbation theory admixtures of small 
components to the chaotic eigenstates, which leads to a non-zero value of S c (|56|), is indeed non-trivial. For example, 
if one examines the sum of Eq. (|56|) as a function of i and j, it would be hard to guess that the sum itself is essentially 

non-zero, since positive and negative values of £jj = J 'cj™ 1 ' 'cj" 2 ' 'cj,™ 2 ^ seem to be equally frequent, and have 
roughly the same magnitude, see Fig. 15. However, in spite of apparent random character of the terms £jj , its mean 



value turns out to be non-zero and is of the same order as the diagonal term Sd ■ Since Sc 



s (mn 2 ) = 

(see below), the suppression of p 2 (a, (3) at the tails should be accompanied by correlational enhancement of the matrix 
elements near the maximum (for — — LOp a \ < T). 

Thus, we come to the important conclusion: even for a random two-body interaction, the correlations produce 
some sort of a 11 correlation resonance" in the distribution of the squared matrix elements p 2 (a,(3). One should note 
that this increase of the correlation effects in the matrix elements of a perturbation can be explained by the increased 
correlations between the Hamiltonian matrix elements when the number of particles and orbitals increases (N/n oc e"). 

In a similar way one can estimate the size of the correlation contribution S c near the maximum of the M 2 distribution 
(at \E^ - -u Pa \< T) ||. 



Sd + Sq 

s d 



1 + Rm =2-(l + R t ) ~2 



n(m — n) 



(62) 



Comparing the values of the ratio S c / Sd at the maximum and at the tail in Fig. 13 (n — 7, m = 14) , one can see that 
indeed, R m rs — Rt- For larger n and to the correlation en hanceme nt factor asymptotically reaches its maximal value 
of 2. This numerical example shows the enhancement of p 2 (a,f3) with respect to Sd at the maximum even greater 
in size than that predicted by Eq. ([32]). This is not too surprising since in Eqs. (pi|)-(|62"|) we estimated the average 
value of R m over an interval AE ~ T around the maximum rather than the peak value at the maximum. 
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FIG. 15. The distribution of the terms = Cf 1 ' 'Cf 1 ' ] Cf 2 ' 'Cf 2 ' of the sum (|56j) for m = 55, n 2 = 66, obtained in the 
TBRI- model for the same set of parameters as in Figs. 1-5, averaged over N g = 100 realizations of the interaction. Indices i 
and j in the figure run over those 84 components in which a is occupied and j3 is vacant, (a) on the top: positive values, (b) 
on the bottom: negative values (absolute values). 



III. THERM ALIZATION AND ONSET OF CHAOS 



A. Distribution of occupation numbers 

Let us now come back to the distribution of occupation numbers defined by Eq.(^|). It gives the probability that 
one of n particles occupies an orbital s specified by the one-particle state \s) , for the fixed exact (compound) state \i) 
. According to this expression, this probability can be found by projecting the state \i) onto the basis of unperturbed 
states, for which the relation between the positions of all particles in the single-particle basis and specific many-particle 
basis state is known by the construction of the latter. One can see that the probability n s = n s (£«) is the sum of 
probabilities over number of basis states which construct the exact state. For Fermi-particles, only one particle can 
occupy an orbital, this is why the occupation number = or 1 . 

It is clear that for chaotic eigenstates the n s — distribution is a fluctuating function of the total energy E — E"> of 

(i) 

a system, due to fluctuations of the components Cf, of eigenstates. In order to obtain a smooth dependence, one 
should make an average over a small energy window centered at E^ , which is in the spirit of conventional statistical 
mechanics for systems in the contact with the thermostat. In fact, such an average is a kind of microcanonical 
averaging since it is done for the fixed total energy E of a system. Therefore, in what follows, by the n s — distribution 
we assume the averaged distribution, 



n s (E) = V Icf 2 (k\n s \k) = y / F(E k ,EW).(k\n s \k) (63) 



k k 

where the F— function discussed in previous Section is used. 

We are going to show that the n s — distribution plays essential role in the statistical approach to finite systems of 
interaction particles. Our interest to this quantity is of two-fold. First, the knowledge of the distribution of occupation 
numbers gives the possibility to calculate mean value of any single particle operator (M) = Yls n sMss ■ Here M ss 
are diagonal matrix elements of a single-particle operator M , which very often can be found easily since they refer 
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single-particle physics. For quantum systems with complex behavior, the non-trivial part is the F— function which 
absorbs the result of the two-body interaction between many basis states. The important point is that in order to 
find the distribution of occupation numbers, there is no need to know exactly eigenstates of the system. Instead, the 
shape of eigenstates in the energy representation is needed, which is defined by the F— function. Thus, if we know 
this function and properties of the unperturbed system, one can relate statistical properties of chaotic systems to 
single-particle quantities. 

Moreover, the variance of the distribution of non-diagonal elements of M , describing transition amplitudes between 
"chaotic" compound states due to a weak external perturbation, can be also expressed through to the occupation 
numbers n s . This variance is important from experimental point of view, for example, for the estimate of an 
enhancement of a weak interaction (which refers, for example, the parity violation in atoms and nuclei) , see details 

Second, the form of the distribution of occupation numbers is interesting itself. One of important questions which 
arises in this respect, is whether the standard Fermi-Dirac and Bose-Einstein distributions occur for isolated systems 
of finite number of particle. If they occur, then, under what conditions and how they can be described? One can 
suggest that the role of the interaction between particles is crucial since this is the only reason to result in a complex 
behavior (chaos) of a system, and the latter is the mechanism for the statistical equilibrium. Another non-trivial 
question relates to the meaning of temperature for isolated systems. The study of the n s — distribution gives a new 
insight on these and other problems. 

For the analytical treatment it is convenient to represent the n s — distribution in the following form, 

J2n<^F(E k -E) 

n s (E) = = (64) 

k 

Here the function F(Ek — E) is the part of the F— function, which non-trivially depends on the difference Ek — E 
between the unperturbed energy Ek only. Namely, we omitted the normalization term / (E) since the summation in 
Eq.(^34j) runs over the unperturbed energy, F{Ek, E) = f (E) F(Ek — E) . Thus, the denominator appears due to the 
normalization. This form (|64| ) allows one to introduce a kind of the partition function, 

Z(E) = Y J F{Ek-E) (65) 
k 

which is entirely determined by the shape of chaotic eigenfunctions. 

The above expression ( |64|) gives a possibility of the statistical description of complex systems. Indeed, as was 
mentioned above, the shape of the F— function has universal features and can be often described analytically. 
Therefore, in practice there is no need to diagonalize huge Hamiltonian matrix of a many-body system in order 
to find statistical averages. We would like again to stress that the summation in (|64| ) is carried out over unperturbed 
energies Ek defined by the mean field, rather than over the energies of exact eigenstates in the standard canonical 
distribution. As a result, the distribution of occupation numbers can be derived analytically (see below) even for few 
interacting particles, in a situation when the standard Fermi-Dirac distribution does not occur. 

B. Microcanonical vs. canonical distribution 

Let us now compare the n s — distribution ( |64| ) with occupation numbers given by the standard canonical distribution 

n, 

^>i 4) exp(-£M/T) 

ns{T) = ^«pW) (66) 

i 

Here T is the temperature of a heat bath and the index i stands for exact eigenstates. The important difference 
between the n s — distribution ( |64| ) and the canonical distribution (|6^) is that in Eq. ( [64] ) the occupation numbers 
are calculated for the fixed total energy E of a system unlike the fixed temperature T in Eq.(|66|). However, both 
expressions can be compared with each other using the relation between the energy E and the temperature T , 



J2E^ exp(-E^/T) 
£exp(-£;tt/T) 



E=(E) T = „ ff . w ^ (67) 
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One can show that the canonical distribution corresponds to the average of the "microcanonicai' n s — distribution 

from Eqs.(§,fl]) into Eq.© 



over some energy interval At . To demonstrate this, let us substitute and 

and replace the summation over i by the integration over p(E^) dE^ where p(E^) is the density of exact energy 
levels, 



C, 



J2 n s ] ex P (~E^/Tj « / n^p(E^) exp (-E® /T\ dE® « 
5>S fc) / Fjp(E k ,E®)p(E®)exp (-E®/t) dE® = J^n?' F(T, E k ) 



(68) 



Here the function F(T,Ek) is the canonical average of Fjp , 

F(T,E k ) = J F^<P T (E^)dE 

with another "canonical (thermal) averaging" function, 

®t(E) = P{E) exp (-E/T) 

As a result, one can transform the canonical distribution ( |66| ) into the form similar to the n s — distribution (|64|), 

,( fe ) F{T,Ek) 



n s {T) 



k 



EF(T,E k ) 

k 



(69) 
(70) 

(71) 



This distribution can be used for the calculation of the occupation numbers and other mean values in a quantum 
dot which is in thermal equilibrium with an environment (with no particle exchange). 

In many-body systems with large number of particles the function $t (E) has a very narrow maximum since the 
density of states p(E^) typically grows very fast. The position E m of its maximum is defined by the expression 



and the width is given by 
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Let us consider the TBRI-model ([[]) for which the density of states is known to be described by the Gaussian, 
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with E c and a 2 as the center and the variance of the spectrum. This allows easily to find the form of <&t(E) which 
also has the gaussian form, 
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One can see that the width Ay of the thermal averaging function is equal to the gaussian width of the spectrum, 
Ay = (t . In Rcf. jl5) it was argued that for large number of particles both widths At and AE are much smaller 
than the typical energy interval, a / \E — E c \ ~ 1/\A^ ■ Therefore, for large number of particles the function $t can 
be regarded as the delta-function at E = E m and the n s — distribution is close to the canonical one, see Eq.(69). 

However, the canonical distribution ( |66| ) is not correct when describing isolated systems with small number of 
particles, instead, one should use the n a — distribution (|64|), see details below and in Jl^Jl5] |. 
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C. Transition to the Fermi-Dirac distribution 



It is naturally to ex pec t that for a very large number of particles the standard Fermi-Dirac distribution arises from 



the n s — distribution (64). Below we reproduce the derivation given in Ref. |15| . By splitting the sum in two parts, 
which corresponds to the separate summation over n s = and n s — 1 , one can represent the expression (|64|) in the 
form 

+ Z a (n-l,E-e a ) _ 1 
s[ >~ ZJn-l,E-i s )+Z s (n,E) i i ^hB 1 ' 

sy ' SJ sv ' ' 1 + Z s {n-1,E-Z B ) 

where two "partial" partition functions Z s (n,E) and Z s (n — 1,E — e s ) are introduced. For the first function the 
summation is taken over all single-particle states of n particles with the orbital s excluded, Z s (n, E) — J^V F(Ek — E). 
Correspondingly, the sum in Z s (n — 1,E — e s ) is taken over the states of n — 1 particles with the orbital s excluded. 
The latter sum results from the terms for which the orbital s is filled (n s = 1) , thus, one should add the energy 
e s = Ek(n) — Ek(n — 1) of this orbital to the energy Ek(n — 1) of the basis state \k) defined by n — 1 particles. Since 
the F— function mainly depends on the difference Ek + e s — E , the adding term e s to Ek(n — 1) is the same as its 
subtraction from the total energy E . Note, that this term is defined by 



Uspn^ (78) 



where e s is the energy of a single-particle state and u sp is the diagonal matrix clement of the two-body residual 
interaction. By taking e s independent of k we assume that the averaging over the basis states near the energy E is 
possible, in fact, this is equivalent to a local (at a given energy) mean field approximation. One should stress that 
this approximation is important when comparing the simple TBRI-model (|l|) with realistic systems. For example, 
for the Ce atom there are orbitals from different open sub-shells (e.g. 4/ and 6s ) which are quite close in energies, 
however, they have very different radius. As a result, the Coulomb interaction between the corresponding electrons is 
very different In this case the interaction terms in Eq.(^|) strongly depend on the occupation numbers of other 
particles. As a result, the equilibrium distribution for occupation numbers n s is very different from the Fermi-Dirac 
distribution p2| . However, the original n s — distribution (|6^ ) for occupation numbers is valid and gives correct result 
p3| . In other cases like random two-body interaction model |T^,|l3|,|l5| or nuclear shell model |2^,^], or the atom of 
gold |54|) , such a local mean field approximation is quite accurate and results in the FD-distribution. 

For large number n 1 of particles distributed over m 3> 1 orbitals, the dependence of Z s on n and e s is very 
strong since the number of terms N in the partition function Z s is exponentially large, N = ^ m ™^ !k! ■ Therefore, to 
make the dependence on arguments smooth, one should consider In Z s instead of Z s . In this case one can obtain 

In Z s (n - An, E - e s ) w In Z s {n, E) -a s An - f3 s e s (79) 
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This leads to the distribution of the Fermi-Dirac type, 

1 



1 + exp(a s + /3 s e s 



(80) 



If the number of substantially occupied orbitals in the definition of Z s is large, the parameters a s and (3 a are not 
sensitive as to which particular orbital s is excluded from the sum and one can assume a s — a = —u/T, (3 S — fj = 1/T 
as in the standard derivation of the Fermi-Dirac distribution for systems in contact with thermostat. Therefore, the 
chemical potential /i and temperature T can be found from the conditions of fixed number of particles and fixed 
energy, 

^n s = n; e s n s + u sp n s n p — n s (e s + e s )/2 = E (81) 

s s s>p s 

Note, that the sums in ( |8l| , |7^) containing the residual interaction u sp can be substantially reduced by a proper choice 
of the mean field basis (for instance, the terms u sp can have different signs in such a basis). In practice, the values e s 
and e s may be very close. Since in the above expressions the nondiagonal matrix elements of the interaction are 
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not taken into account, one can expect that the distribution of occupation numbers defined by these equations gives 
a correct result if the interaction is weak enough (ideal gas approximation). However, below it will be shown that, 
in fact, even for strong interaction the Fermi-Dirac distribution can be also valid if the total energy E is rescaled 
in a proper way, by taking into account the increase of the temperature due to statistical effects of the (random) 
interaction. 

One should also note that the above consideration is similar to the standard derivation (see e.g. J42|) of the Fermi- 
Dirac distribution from the canonical distribution ( |66| ) for the case of many non-interacting particles (ideal gas). It is 
curious that the Fermi-Dirac distribution is very close to the canonical distribution ( |66| ) even for very small number 
of particles, n > 2, provided the number of essentially occupied orbitals is large (which happens for T > e or fi > e ). 
This fact results from the large number of "principal" terms in the partition function Z s , and allows one to replace 
a s by a in the term Z s {n,T)/Z s (n — 1,T) = exp(a s + (3T) in the canonical distribution ( |66| ) (compare with (]77[)). 

One should stress, however, that the temperature T in the Fermi-Dirac distribution is different from that in the 
canonical distribution. Indeed, using the expansion a s — a(ep) + a' {e s — ep) one can obtain the relation between the 
Fermi-Dirac (Pfd) and canonical (/?) inverse temperatures, Pfd = + a'ep ■ Concerning the chemical potential, its 
definition also changes, —fJ,/T = a(ep) — a 1 ■ More specifically, for the same total energy E of the system , the 
canonical and Fermi-Dirac distributions give the same distribution n s defined, however, by different temperatures, 
see details in |I^,|l3|,|l5| and discussion below. 



D. Analytical approach to the n s —distribution 



In the previous section it was shown how the standard Fermi-Dirac distribution occurs in the TBRI-model when 
number of particles is very large. However, the expression (||) for the distribution of occupation numbers via the shape 
of chaotic eigenstates is of more general form and also valid even when the number of particles is relatively small. 
In this case the n s — distribution can be of the form very different from the FD-distribution. Below we show how 
to analytically derive the n s — distribution and express it in terms of single-particle and unperturbed many-particle 
spectrum, using general properties of the F~ function pl| . 

For simplicity, we consider the case of relatively strong interaction, when the shape of the LDOS and exact eigen- 
states can be described by the Gaussian. In order to calculate the occupation numbers n s , we use the expression ( |7^ ) 
containing two partial partition functions Z s (n,E) and Z s (n — X,E — e s ) which correspond to systems with n and 
n — 1 particles, with the orbital s is excluded from the set of single-particle states. The partition function can be 
found from the relation 

Z(E) = F(E k -E)f* J p (E k )F(Ek - E)dE k (82) 

k 

As was discussed above, the density of unperturbed states po(Ek) in the TBRI-model it the Gaussian, s 

where E c is the center of the energy spectrum and N is the total number of states. If the shape of eigenstates is also 
described by the Gaussian, 

F( Ek -E)= - N exp (- {Ek ~ E f \ (84) 
V ; y/MAEY \ 2(A£) 2 J K > 

then the integration in ( |S2| ) can be easily performed. The variance (AE) 2 is defined by Eq. (]3T|) . It should be 
pointed out that, strictly speaking, in this expression the center of the F— function is shifted by the value from 
the unperturbed energy, E — E^> + , see details in p!B[ . This shift is due to the level repulsion which forces 
eigenvalues 2?W in the lower part of the spectrum to move down. The mean-field energies E k = -f/fcfc do not include 
the nondiagonal interaction which results in the repulsion. Therefore, the "center " of the F— function is shifted by 
the value A^ = Ha — E^ . This shift is estimated in Ref. [jl5| as follows, 
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where Uq is the variance of the unperturbed spectrum. One can see that since the variance (AE) 2 of the LDOS is 
typically much smaller tha n ak , this shift in many cases can be neglected. 
Direct integration in Eq.d82) gets 
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(86) 



where a 2 = <Tq + (AE) 2 , therefore, the variance of the partition function Z(E) coincides with the variance of the 
perturbed spectrum. In order to calculate the occupation numbers n Sl one should use the expression @. Therefore, 
one needs to find the partition functions Z s (n, E) and Z s (n — 1,E — e s ) corresponding to n and n — 1 particles, with 
the orbital s excluded from single-particle spectrum. To do this, one needs to calculate the number of states N s and 
the center E cs for these truncated systems, 
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The variance <7q s of the energy distributions can be estimated as 
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°os( n ) ~ CT L n ; <jl s (n - 1) w (cr^) (n - 1) 

where (Tj s is the variance of single-particle spectrum with the excluded orbital s . Here, for simplicity, we have 
neglected the Pauli principle which is valid for to 3> n . Finally, the distribution of occupation numbers takes the 
form 
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where ct? = a 2 Q + (AE) 2 . Numerical data for the TBRI-model are presented in Fig. 16 from which very good agreement 
with (p8|) is seen. 




FIG. 16. Analytical description of the occupation numbers. The data are given for the TBRI-model for the parameters of 
Figs. 1-5 (n — 4, m = 11, Vo = 0.12, do = 1). The histogram is obtained according to (^) by the averaging over eigenstates with 
energies taken from a small energy window centered at E — 17.33 and over 20 Hamiltonian matrices ([j]) with different realization 
of the two-body random interaction. Stars represent the analytical expression ( |8^ ) with ao s found from the single-particle energy 
spectrum. Diamonds correspond to the Fermi-Dirac distribution with renormalized energy, see Section 3.5. 
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It is instructive to compare this result with the Fermi-Dirac distribution which is valid for large number of particles. 



In this case R 
below, see (|9 
particles n 



cxp((e s — /ij/Tth) where T t h = o j(E c — E) is the thermodynamic temperature which is discussed 
The chemical potential /i can be found numerically from the condition of fixed total number of 



E. Effective Fermi-Dirac distribution for finite systems 



In previous Section the distribution of occupation numbers has been derived without any reference to the tempera- 
ture, from the F— function and properties of the unperturbed system. However, the n s — distribution in Fig. 16 seems 
to have a Fermi-Dirac form. One should remind that the latter form in conventional statistical mechanics can be 
derived for ideal gas of very large non-interacting particles. In such a derivation, the presence of the thermostat is 
assumed, actually, in order to have statistical equilibrium in the system. Indeed, for an isolated systems with large 
n — > oo , any extremely weak interaction with an environment results in strong statistical properties of a system. 
Using modern language, one can speak about the onset of chaos due to this interaction. In fact, the weakness of the 
coupling to the heat bath gives the possibility to treat the gas of particles as an ideal gas. It is well known, that in 
this case one can write the following equations, 



E 



n s = n, 



= E 



(89) 



where n and E are total number of particles and total energy, and n s is assumed to have Fermi-Dirac form, 

1 



1 + cxp (a + (3e s ) 



(90) 



When in an isolated system described by the TBRI-model, the number of particles is very large, the above equations 
result in the FD-distribution, see Section 3.3. However, in such a case the interaction Vq has to be very weak. Now, 
if we consider the model with finite and not large number of particles, for a weak interaction there is no chaos in 
the sense that exact eigenstates have small number of principal components N pc -C 1 . Therefore, this model does 
not allow for its statistical description, in other words, there is no statistical equilibrium. For example, if for such a 
case we compute the n s — distribution according to the definition (^), there are very large fluctuations in occupation 
numbers when slightly changing the total energy E^ 1 , see Fig. 17. 




FIG. 17. Absence of a statistical equilibrium for the distribution of the occupation numbers. The histogram is obtained 
in the same way as in Fig. 16, for a very weak interaction V = 0.02 which correspond to the region of (II) of the "initial 
chaotization" , see Section 3.7. The total energy is — 17.33. Stars correspond to the theoretical expression (0) which 
is not valid in this region due to absence of equilibrium. Circles are obtained by direct numerical computation of n s with 
the F— function taken in a specific form, and with the summation performed over real unperturbed spectrum (instead of the 
integration with the Gaussian approximation for p used in Eq. (pq)), see details in |15|. 
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In order to have chaotic eigenstates, and, as a result, the possibility of statistical description, one needs to increase 
the interaction in order to exceed the threshold Vq > df (see also discussion in Section 3.7). The less number of 
particles, the stronger interaction is needed since the two-particle density pf = d^ 1 strongly depends on the number 
of particles. On the other hand, if interaction is strong, the second equation in Eq.(|S9|) is not correct and can not be 
used for the derivation of the n s — distribution. To demonstrate this, in Ref. Jl5| the distribution of occupation numbers 
n s for the two-body random interaction model has been directly computed according to Eq. (|^) from exact eigenstates 
of the Hamiltonian matrix (0), see also H,E^. These data for the "experimental" values of n s are shown in Fig. 18 
by the histogram which is obtained from the average over small energy window in order to smooth the fluctuations 
(with additional averaging over different realizations of the two-body random matrix elements) . To compare with the 
standard Fermi-Dirac distribution, Eqs.(^9|) have been also numerically solved in order to find the temperature and 
chemical potential, the resulting n s — distribution is shown by circles. One should stress that the value of the energy 
E in ( pl| ) was taken the same as for the exact eigenstates from which actual distribution of n s was computed, namely, 
E ~ E^ . The comparison of the actual distribution (histogram) with the "theoretical " one, reveals a big difference 
for a chosen (quite strong) perturbation V — 0.20 . 

To describe correctly the n s -distribution in terms of the Fermi-Dirac distribution, in Ref. |l5[ ] it was suggested to 



renormalize the total energy of the systems due to the interaction between particles, and instead of Eq.(|89[) to solve 
the following equations, 

J2n s = n, ^2e s n s = E + A E (91) 

s s 

where is the shift of the total energy due to the interaction. In Ref. [Jl5| it was argued that in the case of random 
interaction, this term absorbs statistical increase of the energy and gives the correct result for the n s — distribution. 
In fact, this assumption is based on a deep equivalence between the external chaos originated by the heat bath in the 
case of open systems, and internal chaos due to a random character of the interaction. Therefore, it was assumed 
that random interaction and chaos in closed systems plays the role of a heat bath. 




FIG. 18. Fermi-Dirac distribution with and without additional shift of energy due to a (strong) interaction, Vo = 0.20 (other 
parameters are the same as in Fig. 16). Circles stand for the Fermi-Dirac distribution obtained for the total energy E = 17.3 
corresponding to the energy of eigenstates, see (|si|). Diamonds correspond to the distribution obtained for the energy shifted 
according to Eq. (p^). 



In order to find analytically the shift A e , one needs to consider the structure of exact eigenstates in the unperturbed 
basis, in particular, to find the shift between the energy of an exact eigenstates and the mean energy of the components 
of the same eigenstate []l5| . Since the density of states rapidly increases with the energy E , the number of higher 
basis states admixed to an eigenstate due to the interaction is larger than the number of lower basis states. As a 
result, the mean energy 
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(E k ) i = J2 E k F^ « f E k F®po(E k )dE k (92) 

k 

of the components in an exact eigenstate \i) is higher than the eigenvalue E^ 1 corresponding to this eigenstate (we 
consider here the eigenstates in the lower part of the spectrum). There is another effect which decreases the value of 
(Ek) i , see (|92l), which remains even if the density of states does not depend on the energy. This second effect is due to 
repulsion between energy levels, according to which the eigenvalues move down for this part of the spectrum, therefore, 
the difference between {Ek) i and E^' increases due to the interaction. The second effect also shifts the "center" of 
the F— function. One should stress that all effects leading to the above shift of the energy are automatically taken 
into account in the relation (|92]). Thus, one can analytically calculate this shift A# = (Ek) i — E^ from the equation 
(|92|). For this, one needs to know the unperturbed density of states and the form of F— function. The evaluation of 
the shift Ae has been done in fig] by assuming some generic form for the F— function which is valid in a wide range 
of the interaction strength V , 

A E = (E k ) l - E« = ^£ (E c - E) (93) 
a o 

where E c is the center of the energy spectrum. 

Thus, to find correct values for the occupation numbers in the Fermi-Dirac distribution, we should solve Eqs. 



with defined by Eq.(93). The resulting n s — dependence is shown in Fig. 18 by diamonds. As one can see, sue 
a correction gives quite good correspondence to the numerical data. Similar difference occurs for larger number of 
particle (and smaller interaction strength), see Fig. 19 where due to serious numerical problems, the data for the 
occupation number distribution are given without direct comparison with actual n s — distribution. 




FIG. 19. The same as in Fig.18 for large number of particles, n = 14, and orbitals, m = 28 , for Vo = 3 • 10 3 and do = 1 
. The histogram is not possible to obtain numerically in the direct computation, due to a very large size of the Hamiltonian 
matrix, N = 330000 . 



To check the analytical prediction (£3|) for the shift A^ , this shift has been directly calculated by comparing 

the energy E^ of exact eigenstates with the energy (E k ) i . The latter has been numerically found from the exact 

2 



relation (Eu 



Efc^fc 



(compare with (|92|)). The comparison of these data (circles in Fig. 20) with Eq. 

(straight full line) shows reasonable agreement, if to neglect strong fluctuations around the global dependence. These 
fluctuations are due to fluctuations in the components of specific exact eigenstates \i) (note, that the presented data 
correspond to individual eigenstates, without any additional averaging). 
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FIG. 20. Shift of the total energy for the corrected Fermi-Dirac distribution. The data are given for the TBRI-model with 
n = 4, m = 11, do = 1, V = 0.12 . The straight line is the analytical expression (p3|), the dotted line (squares) presents direct 
computation of the shift based on the diagonalization of the Hamiltonian ([l]) with the following computation of the < Ek >i ■ 
On the horizontal axes the rescaled energy = {E^ — Ef ermi )/(E C — Ef ermi ) is plotted. 



F. Temperature vs. chaos 



In previous Section it was shown that if the n s — distribution for finite number of particles in the TBRI-model has 
the from of the Fermi-Dirac distribution, it can be found from the modified Eq. (^Tj) . These equations also determine 
the parameters a and T = /9 _1 which may be associated with the "chemical potential " and "temperature " for an 
isolated system. One can see, that due to the shift of the total energy in Eq.(|9l|), there is a corresponding shift 
of the temperature AT . It is important to note that in systems with infinite number of particles, all definitions of 
temperature give the same result. In contrast to this, for finite number of particles different definitions give different 
results. 

Let us, first, start with the standard definition of temperature, 

1 = dS th = din p 
T th dE dE 1 ' 

where Sth is the thermodynamical entropy, 

S th = In p{E) + const (95) 

and p{E) is the density of states. 

In fact, this definition of the thermodynamical temperature stems from the estimate of the position of maximum of 
the canonical averaging function <&t (E) , see ([70|) , if we assume that the position of its maximum E m coincides with 
the energy E of a system. One should stress that in the above definition p(E) is the total density of states, therefore, 
the interaction is essentially taken into account. 

Another definition, which is consistent with the first law of thermodynamics (energy conservation), is given by the 
relation (E) T = E , see Eq.(|67|). Here the averaging is performed over the canonical distribution (|66|). Since the 
width At of the canonical averaging function &t{E) is not zero, the two definitions of the temperature, (|94D and ( |67| ) 
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give, in principal, different values. Indeed, in the case of the Gaussian form of p(E) the value of T t h given by ( |94| ) 
takes the form (see also |^6| , ^8| ), 

T - = E^E (96) 

where E c and a are the center and the width of the total density p(E). 

On the other hand, direct evaluation of the relation (pfl) leads to the following definition of the canonical temperature, 



a 2 



E c - E + A 

Here, the shift A is approximately given by the expression 



- = T< 
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(97) 



(98) 



where 



K= I rx P ( ~) dx « v^F ; x = - Em ; E m =E C -^— (99) 

^ / ^ ^ can 

One can see that the shift A itself depends on the temperature, it is proportional to the width At =o~ of the function 
<$>t{E). In the above expressions, -Emin and E mSuX are the low and upper edges of the energy spectrum. Note that 
the relation A = occurs at the center of the spectrum, therefore, at the center the temperature diverges and in the 
upper part of the spectrum is negative. This fact is typical for systems with bounded spectrum, for example, for spin 
systems. In fact, the TBRI-model (0) with finite number m of orbitals can be treated as a model for one open shell in 
atoms, nuclei, clusters, etc. However, in realistic many-body systems there are always higher shells which contribute 
to the density of states for higher energy. Thus, the density of states p(E) is a monotonic function which results in 
the positive temperature. For such physical applications, the model (|l|) with finite number of orbitals is reasonable 
in the lower part of the energy spectrum where the influence of higher shells can be neglected. 
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FIG. 21. Different temperatures versus the rescaled energy x = (E—Ef erm i)/(E C —Ef erm i) for the TBRI-model . Parameters 
are the same as in Fig.16. Triangles stand for the thermodynamical temperature Tth denned by ( |96| ) and should be compared 
to the canonical temperature T can (circles), see ([^)- The width a of the perturbed density of states is defined by the residual 
interaction Vb = 0.12 according to (^l|) and the relation a 2 = a\ + (AE) 2 with <jo found numerically from the unperturbed 
many-particle energy spectrum. 
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One can also see that the difference between the two equations of state T{E) defined by Eqs.(|96|) and (J97|), disappears 
for highly excited eigenstates (for which E m — E m i n 3> cr), or in large systems with n 1 . Indeed, one can obtain, 
E c — E ~ no\ , where o\ is the width of single-particle spectrum. On the other hand, according to the central 
limit theorem, the variance of total energy spectrum can be estimated as dp = ^ n a\ sa naf , therefore, the ratio 
a/(E c — E) ~ 1/s/ri tends to zero at n — ► oo . As was mentioned above, in such realistic finite systems like atoms 
and nuclei, the number of particles in an open shell is relatively small (n — 4 for the Ce atom |19| and n = 12 in 
nuclear shell model ]2^ , |28[ |), therefore, the corrections to the thermodynamical temperature ( p4| ) may be important, 
especially, for low energies. The detailed discussion of different temperatures in nuclear shell model is given in |p6| , p8| . 

The energy dependence of temperatures T t h and T can , as well as the temperature T exp found directly from the 
numerical simulation, is shown in Fig. 21. The data refer the TBRI-model with n — 4 interacting Fermi-particles and 
77i = 11 orbitals. The comparison of the thermodynamical temperature T t h defined by (|96}), with the "canonical" 
temperature ( |97| ) reveals quite strong difference in all the range of the rescaled energy \ = (E — E /ermi) / (E c — E fermi) ■ 

Now let us find the shift of the temperature AT which is due to the interaction, see previous Section. Since it is 
directly related to the shift of total energy, E = {Ek) i = E™ + Ae , one can get, 

T = T + AT = 3° a ^ ... + {AE) ]., (100) 

Therefore, statistical effects of random interaction can be directly related to the increase of temperature of a system, 
AT/T = {AEf/al . 

One of the important questions deals with thermodynamical description of isolated systems of interacting particles. 
In any thermodynamical approach one needs to define, in a consistent way, such quantities as entropy, temperature 
and equation of state. This problem has been recently discussed p6| , p8| in application to shell models of heavy nuclei. 
In particular, it was shown that for a realistic residual interaction different definitions of temperatures give the same 
result. 



G. Transition to chaos and statistical equilibrium 

Let us now summarize the discussed above results for the TBRI-model in what concerns transition to chaos and 
thermalization. The latter term is not defined for isolated systems, our suggestion is to treat "thermalization " as the 
onset of statistical equilibrium. The latter allows to give a reasonable statistical description and may be used to find 
thermodynamical description. 

Depending on the strength of (two-body) interaction between particles , one can fix different situations in the 
model. First region (I) refers strong (perturbative) localization. This occurs when the interaction is very weak, 
Vq <C df and standard perturbation theory gives correct result. In this case exact eigenstates have only few relatively 
large components (N pc ~ 1) , in other words, the eigenstates are strongly localized in the unperturbed basis. This 
situation is quite typical for lowest eigenstates (where the density of states is small) even if for higher energies the 
eigenstates can be considered as very "chaotic" ones (N pc 3> 1) . 

The second region (II) is characterized by an initial chaotization of exact eigenstates which corresponds to a 
relatively large, N pc 1 number of principal components and V < df. The latter condition is essential since it 
results in very strong (non-Gaussian) fluctuations of components Cjp [^5| for the fixed energy i?W of compound state 
\i) . Such a type of fluctuations reflects itself in a specific character of eigenstates, namely, they turn out to be sparsed. 
As a result, the number of principal components can not be estimates as N pc ps T/D , as is typically assumed in 
the literature. In this case the energy width T of both the LDOS and eigenstates can be approximately describe as 
r « 2irV 2 /df . In Fig. 17 it is shown how the distribution of occupation numbers n s looks like for the TBRI-model (jl]) 
with n — 4 particles and m = 11 orbitals and very weak perturbation V/do ps 0.02 . One can see that the distribution 
of occupation numbers has nothing to do with the Fermi-Dirac distribution (full diamonds), it turns out to be even 
the non-monotonic function of the energy e s of orbitals (see also [^2|). Note that the averaging procedure used in 
Fig. 17 can not wash out strong fluctuations in occupation numbers n s . 

With further increase of the interaction, where N pc 3> 1 and V > df , the region (III) of the statistical equilibrium 

emerges. In this region the fluctuations of eigenstate components are of the Gaussian form |^5| and one can 
introduce the F— function ( |64| ) as the shape of exact eigenstates in the unperturbed energy basis. Correspondingly, 
the fluctuations of the occupations numbers n s are small in accordance with the central limit theorem, An s /n s ~ 

— 1/2 

N pc <C 1 for n 8 *~ 1 . One should stress that in this region the value of N pc is given by the common estimate, 
N pc ~ T/D . As a result, the n s — distribution changes slightly when changing the energy of a system. Such a 
situation can be naturally related to the onset of thermal equilibrium, though the form of the distribution n s can be 
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quite different from the Fermi-Dirac distribution. In this case, the F— distribution allows for a correct description of 
an actual distribution of occupation numbers in isolated quantum systems of interacting particles. One can see that 
the equilibrium distribution for the occupation numbers arises for much weaker condition compared to that needed 
for the Fermi-Dirac distribution. Since the energy interval df between directly coupled basis states is small, it is 
enough to have a relatively weak residual interaction V > df in order to have the equilibrium distribution (note, that 
the value of df decreases rapidly with the excitation energy) . 

Next region (IV) is where the canonical distribution ( pq ) occurs; for this case in addition to the equilibrium, one 
needs to have large number of particles, n 1. If, also, the condition T <C ndo is fulfilled, the standard Fermi-Dirac 
distribution is valid with a proper shift of the total energy due to the interactior^see Section 3.5. Typically, this region 
is associated with the onset of the canonical thermalization (see, for example, 26,281]). In practice, the condition (IV) 
of the canonical thermalization is not easy to satisfy in realistic systems like atoms or nuclei since n in the above 
estimates is, in fact, the number of "active" particles (number of particles in a valence shell) rather than the total 
number of particles. Thus, the description based on the F— distribution (64) which does not require the canonical 
thermalization condition (IV), is more accurate. 

The above statements are confirmed by the direct numerical study of the two-body random interaction model []l2|jl5| 
with few particles when changing the interaction strength V/dg. If, instead, we increase the number of particles keeping 
the interaction small, V <C do , the distribution (|64|) tends to the Fermi-Dirac one as it is expected for the ideal gas, 
see El. 



IV. CONCLUDING REMARKS 



In this paper we have discussed a novel approach to isolated systems of finite number of interacting particles. 
The goal of this approach is a direct relation between the average shape of exact eigenstates (F— function), and 
the distribution n s of occupation numbers of single-particle levels. From this relation one can see that there is no 
need to know exactly the eigenstates, instead, if these eigenstates are chaotic ( random superposition of a very large 
number of components of basis states), the F— function absorbs statistical effects of interaction between particles and 
determines the form of the n s — distribution. Therefore, the structure of chaotic eigenstates in dependence on the 
model parameters, is the central question in this approach. 

The results discussed above relate to the TBRI-model for which all two-body matrix elements are assumed to be 
random and independent variables. This assumption was made in order to study limiting statistical properties of 
this model. In particular, it was shown that even in this limit case of completely random (two-body) interaction, 
the Hamiltonian matrix in many-particle representation can not be treated as the random matrix, therefore, the 
RMT is, strictly speaking, not valid. However, under some conditions exact eigenstates turn out to be quite random 
and statistical approach is correct, however, one should take into account the form of eigenstates in a given basis of 
unperturbed part Hq . 

Another reason for the study of this TBRI-model, is that it is very convenient for the demonstration of the 
developed approach. In this case there are no any effects of regular motion in the system, and many of analytical 
estimates can be obtained in a clear way. In particular, it was shown how to calculate the n s — distribution from the 
shape of eigenstates, provided the unperturbed spectrum of energy is known. One can stress that in this way one 
can analytically obtain the n s — distribution which can have the form very different from the standard Fermi-Dirac 
distribution. On the other hand, for sufficiently large interaction the n s — distribution is of the Fermi-Dirac form, 
therefore, it is convenient to introduce an effective "temperature " and "chemical potential " which give correct 
description of actual n s — distribution in terms of the Fermi-Dirac distribution. To do this, one needs to find the shift 
of the total energy which formally comes into equations determining the FD-distribution. One should stress that this 
shift is directly related to the F— function and can be found analytically. 

Now, we would like to point out that the approach discussed in this paper for the TBRI-model is of generic and 
can be applied both for random and dynamical interaction. One of the most interesting problems is the application 
of the approach to dynamical systems with the well-defined classical limit. In this situation exact eigenstates in the 
corresponding quantum model can be expected to appear when two conditions are fulfilled. The first one is the strong 
chaos in the classical counterpart, and the second is the semi-classical limit (which is typically equivalent to a high 
energy of a system). Under these two conditions, eigenstates of quantum model have many components and these 
components may be treated as pseudo-random, thus leading to a statistical equilibrium in the system and possibility 
to apply the suggested approach. 

In Refs. |^,[l8|,[l7| two quantum dynamical systems have been studied and compared to their classical limits. In 
both cases the main results refer the region of parameters where the classical motion is strongly chaotic. One of the 
important questions which was under close investigation is the quantum-classical correspondence for the F— function 
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and the LDOS. As was pointed out in Ref. |39[] , there is a quite clear and easy way for finding classical F— function 
and classical LDOS. It is instructive to explain this approach since it is of generic and can be used in many physical 
applications (see details in |40UT^ , |l7t ) . 

Let us start with the classical F— function. We assume that the total Hamiltonian can be represented in the form, 

n 

H = H + V ; H = J2 H k; H k ° = H°( Pk ,x k ) (101) 
fc=i 

Here Hq stands for the "unperturbed " Hamiltonian which is the sum of partial Hamiltonians H k describing the motion 
of different (non-interacting) n particles. The interaction between particles is embedded in V which is assumed to 
result in chaotic behavior of the (total) system. Note that the same consideration is valid if instead of particles we 
mean different degrees of freedom for one particle. Now let us fix the total energy E of the Hamiltonian H(t) and 
find (numerically) the trajectory Pk{t) , Xk(t) by computing Hamiltonian equations. Since the total Hamiltonian is 
chaotic, there is no problem with the choice of initial conditions Pfc(0) , 2^.(0) , any choice gives the same result if one 
computes for sufficiently large time. When time is running, let us collect the values of unperturbed Hamiltonian Ho(t) 
for fixed values t — T, 2T, 3T, ... , and construct the distribution of energies E (t) along the (chaotic) trajectory of 
the total Hamiltonian H . In such a way, one can get some distribution W(Eq;E = const) . Comparing with the 
quantum model, one can see that this function W{Eq\E — const) is the classical analog of the F— function which 
is the average shape of eigenstates in energy representation. Indeed , any of exact eigenstates corresponds to a fixed 
total energy E = const and it is represented in the unperturbed basis of Ho , in fact, F~ function is the (average) 
projection of exact eigenstate onto the set of unperturbed ones. Thus, one can expect that for chaotic eigenstates in 
a deep semiclassical region the two above quantities, classical and quantum ones, correspond to each other. 

On the other hand, one can consider the complimentary situation. Let us fix the unperturbed energy E and compute 
a trajectory p k °^ (t) , x k (t) which belongs to the unperturbed Hamiltonian Ho(t) . Similar to the previous case, let us 
put this unperturbed trajectory into the total Hamiltonian Hit) and collect the values of total energy E(t) along the 
unperturbed trajectory for discrete values of time. In this way one can find the distribution W{E; Eq = const) which 
now should be compared with the LDOS in the corresponding quantum model. However, in this case one should be 
careful and make an average over many initial conditions Pfc(0) , x k (0) with the same energy Eo , if the unperturbed 
Hamiltonian is regular. In fact, the above two classical distributions W(Eq\E = const) and W(E;Eq = const) 
determine the ergodic measure of energy shells, the first one, when projecting the phase space surface of H onto H a 
, and in the second one, the surface of Ho onto H (see discussion in ]39[). 

Numerical data for two different dynamical models have shown amazingly good correspondence between the 
F— function and the LDOS, and their classical counterparts, see details in Refs. p0| , |l8| , |l7| | . Recently, few other 
systems have been studied both in classical and quantum representations, and again, very good correspondence has 
been numerically found in the semiclassical region. These data confirm the theoretical predictions, and seems to 
be very important in view of future developments of the semiclassical theory of dynamical chaotic systems. It is 
important to point out that the above correspondence opens a new way for the semi-quantal description of quantum 
system, when the form of the F— function is taken from a classical model and used in order to find distribution of 
occupation numbers n s of single-particle levels in the corresponding quantum system |56[ | . 

Another interesting problem which has been studied in dynamical systems, is the distribution of occupation numbers 
and the possibility of its analytical description in the same way as it was done for the TBRI-model (see details in 
[p7t). One of the most interesting results obtained numerically, is that the canonical distribution occurs in an isolated 
(dynamical) system of only two interacting spin-particles, if one randomize the non-zero elements of the interaction 
V as close as possible to the dynamical constrains of the model. This means that random interaction indeed plays a 
role of the heat bath and allows to use statistical and thermodynamical description for isolated systems. 
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